A Fixed-Point Iterative Method for Discrete Tomography Reconstruction Based on Intelligent Optimization

Luyao Yang,Hao Chen,Haocheng Yu, Jin Qiu,★and Shuxian Zhu,★

1School of Electronic and Information Engineering,Suzhou University of Science and Technology,Suzhou,215009,China

2College of Mechanical Engineering,Suzhou University of Science and Technology,Suzhou,215009,China

ABSTRACT Discrete Tomography (DT) is a technology that uses image projection to reconstruct images.Its reconstruction problem,especially the binary image(0–1 matrix)has attracted strong attention.In this study,a fixed point iterative method of integer programming based on intelligent optimization is proposed to optimize the reconstructed model.The solution process can be divided into two procedures.First,the DT problem is reformulated into a polyhedron judgment problem based on lattice basis reduction.Second,the fixed-point iterative method of Dang and Ye is used to judge whether an integer point exists in the polyhedron of the previous program.All the programs involved in this study are written in MATLAB.The final experimental data show that this method is obviously better than the branch and bound method in terms of computational efficiency,especially in the case of high dimension.The branch and bound method requires more branch operations and takes a long time.It also needs to store a large number of leaf node boundaries and the corresponding consumption matrix,which occupies a large memory space.

KEYWORDS Discrete tomography;integer programming;fixed-point iterative algorithm;intelligent optimization;lattice basis reduction

1 Introduction

The purpose of tomography is usually to reconstruct some physical objects through projection.The reconstruction algorithm plays an important role in the fields of medicine,image processing and computed tomography.Discrete tomography(DT)is essentially a function formed by reconstructing discrete values from multiple projection directions [1].DT regards the reconstructed object as a limited set of pixels, and the slope of the straight line along the projection direction should be a rational number.The framework studied in this work is a two-dimensional lattice, that is, the reconstructed object is an image.From a long history, this field is gradually developed by multiple branches of mathematics.The combination of binary matrices is one example.Ensuring that the reconstructed result is equal to the original object is an important task of tomographic reconstruction.A core problem of DT is to reconstruct (0–1) discrete matrix from multiple directions Herman et al.[2] introduced relevant applications and theories for this problem.Weber et al.[3] focused on reconstructing binary functions, this approach extends the linear programming relaxation of combinatorial optimization problems proposed by previous researchers from a small number of X-ray projections from a small number of X-ray projections to the objective function of quadratic smootha priori,and can calculate the solution that biases the reconstruction to the region with spatial coherence by any interior point solution method.The image reconstruction problem forms a 0–1 discrete matrix through ray projection,which can be expressed in the form of integer programming.

Linear integer programming model is widely used in combinatorial optimization and emergency management.Haghani et al.[4]modeled the emergency vehicle scheduling as an integer programming model to solve the problem of reasonable vehicle routing in the case of serious traffic congestion,and the vehicle routing was flexibly scheduled using real-time traffic data by studying the optimization model.Considering that the rescue route is blocked due to various uncertain traffic accidents, Chai et al.[5] and others proposed a reliable integer programming model to calculate the travel time required for the whole rescue route and optimize the allocation of emergency resources.Cao et al.[6]designed an emergency organization optimization strategy based on branch and bound method to solve the sustainability problem in the disaster response decision optimization system.Huibregtse et al.[7]designed and developed an optimized evacuation instruction framework by using a fixed-point formula, to better optimize the evacuation efficiency of traffic.Dulebenets et al.[8] constructed an integer programming model to guarantee smooth evacuation of people and avoid congestion when disasters occur,and it was combined with heuristic algorithm to reduce the total traveling salesman time of individual evacuation to a certain extent,for optimizing the final evacuation plan.French et al.[9]avoided floating-point data calculation by using fixed-point algorithm,and only processed integer data,to reduce the processing time in aerial positioning and optimize the imaging system.

This reconstruction problem is a NP hard problem.Graph partition[10]and quadratic knapsack problem[11]are also classical NP combinatorial optimization problems.Dang et al.[12,13]proposed a fixed-point iterative method to solve such problems.This method can find an integer point in any type of polyhedron.This fixed-point algorithm has been widely used in the reassignment of aircraft disturbed by aviation.A large number of flights are disturbed and cannot operate according to the original plan due to the emergency bad weather.Wu et al.[14] and others used the distributed integer programming to re-optimize the aircraft allocation based on the fixed-point iterative method for calculating the feasible route.With a slight modification of this method, all integer points in the polyhedron can be further obtained, and then multiple feasible solutions can be obtained.At present, this method has effectively solved many integer programming problems.When solving large-scale integer programming problems, this method can better reflect its effectiveness compared with the previous traditional methods.Our future work can focus on the implementation of fixed-point iterative method with distributed computing,and distributed problems have been widely used in some aspects,such as distributed permutation flow shop scheduling problem.A multi-point iterative greedy algorithm is proposed in [15] to solve this problem.If the problem is not improved after a certain number of continuous iterations, then the restart strategy is adopted to minimize the factory completion time.A mixed integer linear model is established in [16], and the search ability is improved by an improved artificial bee colony algorithm.The heuristic method based on greedy rules proposed in[17]can also effectively solve this problem.Therefore,distributed computing will occupy an important position in the future work.

In the process of solving the DT problem, if its reconstruction problem can be rewritten into the form of integer programming, this problem can be effectively solved by the fixed-point iterative method.On the basis of algorithm proposed by Aardal et al.[18,19],we transform the DT problem into an integer programming system.Branch and bound method is a classical method, which can solve integer programming systems.Its characteristic is that it can quickly obtain the optimal solution.However,the numerical results show that the fixed-point iterative method is more effective and feasible than the two methods mentioned above especially in the case of high dimension.

The fixed-point algorithm proposed in this study has been widely used in the rapid response of emergency management in smart city,such as flight delay,fire and other emergencies.On the basis of this fixed-point theory, an integer programming algorithm is proposed to quickly solve the integer programming of new scheduling under emergencies for providing support for the development of smart city and its core algorithms.It can also provide a new idea and method for the research of integer programming.

The rest of the paper can be divided into the following sections.Section 2 introduces the problem and model of DT reconstruction.In Section 3, LLL(Lenstra-Lenstra-LovaszLattice) algorithm is briefly described, and DT problem is transformed into polyhedron judgment problem through this algorithm.Section 4 shows the calculation principle and workflow of Dang fixed point iterative method.In Section 5,the experimental data are obtained by running MATLAB and then compared with those of the branch and bound method.Section 6 summarizes this work and puts forward the future prospects.

2 Problem Statement and Model

DT essentially re-establishes the function with known discrete range, which is based on the projection along a certain direction.Reconstruction of discrete(0–1)matrix from multiple projections is a key step in the study of DT problem[20].

The following model is more intuitive and clearly shows the essence of the reconstructed object.This model does not involve any prior knowledge related to the reconstructed object.In general,the reconstruction results produced by this model are not the only solution.

Letf:Ω ⊂R2→R,yf (y)be an unknown attenuation function.This study establishes a discrete model similar to that in Fig.1, which can more intuitively describe the reconstruction off.The functionfcan be represented by its valuexj∈{0,1},j= 1,2,···,n,x=(x1,x2,···,xn)T, in a square pixel.We definef (y)=∑xjΦj(y)= 〈x,Φ(y)〉, Φ(y)=(Φ1(y),Φ2(y),···,Φn(y))T,where if and only ifyis within thejth pixel range,then the value of Φj(y)is 1 and 0 at other positions.Several raysi= 1,2,···,mintersect with Ω, the contributionaijof each pixeljforms a projection matrixAofm×n,and theith componentbiof vectorbrepresents the total attenuation along theith ray.In this study we ignore any“noise”generated in the imaging process,and we default the projection ray,point source and detector to the ideal state.

The imaging process is shown in Fig.1.The DT problem can usually be simplified to a 0/1 integer linear programming system in the following form:

We findx∈{0,1}nthat satisfies Eq.(1),where the coefficient matrixA∈Nm×nof the equations is composed ofmrows andncolumns, andb∈Nmis an m-dimensional column vector.The model described above is widely used in the research of material science.For example,when studying crystal materials, the lattice formed by atoms of crystal solids can be transformed into integer point lattice through modeling.Then, obtaining all atoms on a line by using electron microscope technology is relatively easy, and the calculated number can be used as the component of vectorbto participate in the operation of (1).Therefore, after the abovementioned integer linear equations are solved, the problem of image reconstruction can be solved.

Figure 1:Discrete model of DT reconstruction problem[3]

However, in the actual calculation process, the solution of (1) is difficult to find smoothly.The scale of this kind of problem is NP hard.In fact,nearly no suitable solution to solve this problem is available in practical application.

The problem of n-dimensional scale corresponds to the number of atoms contained in the reconstructed object,which is a huge number that is difficult to calculate,regardless of scene application.For example,the most common daily water contains at least 1023 atoms per milliliter.Previous research and calculation show that relevant algorithms are available for solving integer linear programming,but they are time consuming and costly.

Multiple solutions of(1)may be available,but only one solution matching is possible in the actual reconstruction object.Usually,some optimization algorithms are used to solve this problem,but this way increases the computational difficulty of the problem to a certain extent.The solution obtained by calculation can only approximate the actual value in any case.A more accurate approximation can be achieved by increasing the number of projection directions,but this way will undoubtedly increase the scale of the problem again.

This study presents a new solution to problem(1),which avoids the abovementioned difficulties to a certain extent.The purpose is to clarify the computational theory in DT.The characteristics and basic idea of the new algorithm(based on 2D)are as follows:

1.Any possible situation is included in the established calculation model.

2.The projection ray only occurs in one direction with a certain slope,such that only one integer point can be generated on this line at most.

3.Each projection ray to be solved is expressed as linear Diophantine equationax+by=c,wherex,y∈Z,a,b,c∈R.The coefficientsaandbin different linear equations are proportional to each other or the same,that is,the slope of the straight line is the same.

All the solutions obtained by calculation are the accurate reconstruction of the processed object.

3 Transformation of the DT Problem into a Polyhedron Judgement Problem

In [18,19], an LLL [21] based reduction method is proposed, which can transform problem (1)into a polyhedron judgment problem.

If a basec1,c2,···,chofRhexists in subsetLof n-dimensional vector spaceRh,wherehis a positive integer,then it can be called lattice,such that

In this way,c1,c2,···,chcan be called a base ofL,and the rank ofLish.Then the determinant ofLcan be defined as the following form:

wherecjis expressed as a column vector.The result of this determinant is a positive real number,which is independent of the selection of the base.

Given that the relationship betweenc1,c2,···,chis linearly independent,these linearly independent column vectorscj,1 ≤j≤hare orthogonalized into,1 ≤j≤hby an algorithm for deriving orthogonal vectors,that is,Gram-Schmidt orthogonalization rule.MatricesC=(c1,c2,···,ch)andare composed of column vectorscjandc*j,respectively.The definition of vectorcan be summarized as follows:

where

where(,)represents the inner product of two ordinary vectors onRh.The relationship betweenc*jandcjis that the former is the orthogonal projection of the latter onj≤h.Therefore,through the operation of the abovementioned rules,the obtainedc1,c2,···,chis the orthogonal basis onRh.

Notably,C*is not the base ofLin general, because all the values ofμkjare not integral.The following is the definition of approximate orthogonal basis.

Let a set of linearly independent vector groupsC=(c1,c2,···,ch)be a basis ofL,and a new set of basesis obtained through Schmidt orthogonalization.If the following conditions are met,thenCis called reduced basis.

and

‖‖in the formula represents the ordinary Euclidean length.The constantin the inequality is not necessarily a fixed value,it can be randomly selected,and its value can be replaced by any other real number

To obtain the reduced basis of the lattice given according to the original basis,a polynomial time algorithm is introduced below.Its components are size reductions and interchanges.

Size reduction:If the condition of inequality(6)is not satisfied for any indexkandj,1 ≤k≤j≤h,the value ofcjcan be re expressed bycj-ck,whereis the approximate integer rounded byμjk.

Interchange:Ifjfor some indexes does not meet the requirements in(7),the vectorscjandcj+1can be swapped.

After the abovementioned description,the process of base reduction algorithm can be summarized as follows(Fig.2).

Figure 2:Base reduction algorithm

Based on the abovementioned base reduction algorithm, given two suitable and relatively independent positive integersN1andN2,they should be sufficiently large compared with the input.

Then we consider the following matrix,in which any column vectors are linearly independent:

whereI(n)represents n-dimensional identity matrix, 0(n×1)represents an zero matrix, which isn×1-dimensional,Arepresentsaijin the reconstructed model, andbrepresentsbiin the abovementioned model.For latticeB, a new matrix is obtained by applying the LLL basis reduction algorithm of polynomial,which can be re expressed as the following formula:

If a vector exists,then

Through the above analysis,we successfully transform problem(1)into problem(10),that is,to judge whether there are integer points exist in the polyhedron is determined.We use the fixed-point iteration method[9,10]to solve it.If integer points exist,then a solution exists;otherwise no solution exists.We will describe this method in detail in the next section.

4 The Fixed-Point Iterative Method

Let

whereD∈Ra×bis ana×bmatrix withb≥2,where the elements are integers,G∈Ra×pana×pmatrix,andda vector ofRa.In this study,we assume thatPis a region with wide dimensions and upper and lower bounds.For a real numberα,we useto represent the largest integer rounded to the left ofα.Forx=(x1,x2,...,xb)T∈Rb,let

Letxmax=(x1max,x2max,...,xbmax)Twithxjmax= maxx∈P xj,j= 1,2,...,b, andxmin=Thenxmin≤x≤xmaxfor allx∈P.Letwhere,therefore,all integer points,which are inP,are inD(P).In the general case,we assume thatxl ＜xmin(letxil=ximin-1 ifxil=ximinfor somei∈N).Forz∈Rbandk∈N0,letp(z,k)={x∈P|xi=zi,1 ≤i≤k and xi≤zi,k+1 ≤i≤b}.

An integer point is given as the initial pointy∈D(P)withy1＞x1l;We introduce a method to judge whether an integer pointx*∈Pwithx*≤l yexists,which is the fixed-point iterative.The pseudo code for a particular algorithm is illustrated in Fig.3.

Figure 3:Pseudo code of specific algorithm

The basic idea of solving the abovementioned integer problem involves the incremental mapping,which is from finite lattice to itself.Except for the integer points in polyhedronP, all other integer points can be mapped to the first point inP,which is smaller than them in the lexicographical order ofxl.After the abovementioned incremental mapping steps,all integer points in the polyhedron are fixed points.We present an integer point as the initial point,find an integer point in the polyhedron by this method,or prove that no such point exists in the polyhedron in finite iterations.

5 Numerical Results

After the analysis in the second section, we can conclude that the DT problem is a special case of the market split problem[22].The problem of market segmentation is defined as that a company supplies several products to retailers,and allocates retailers to the two departmentsD1andD2in the company.ThusD1can grasp 40%of the market occupied by each product of the company,andD2can grasp the remaining 60%.If the product cannot be divided in proportion completely and accurately,then the proportion error of the division will be minimized.The market segmentation problem can be expressed by a set of mathematical methods of 0–1 linear programming:

In this section,we take market split problem as an example to show the relevant numerical results of solving problem(1).Given any real numberm,the expression ofnisn=10(m-1),the elements in the coefficient matrixAare integers randomly selected from[0,99],and each component of vectorbis calculated by1 ≤i≤m.LetN1=1000 andN2=10000 in the computation.Through the abovementioned examples,we compare the method introduced in this study with the branch and bound method, which are two classical algorithms to solve integer programming problems.We also analyze the number of iterations and time required by the two algorithms to solve the problem.

The design idea of the fixed-point algorithm is to give a polyhedron and define an incremental mapping from the problem lattice to the lattice itself in combination with the incremental mapping principle.In that way, the integer points outside the polyhedron are mapped to the feasible integer points in the polyhedron that are smaller than the dictionary order of the integer points.When the feasible solution is not satisfied,the next round of iteration is conducted.The algorithm is terminated until a feasible solution is found or no feasible solution is found.The basic idea of branch and bound method is to select different branch variables and sub-problems for branching.Usually, all feasible solution spaces are gradually decomposed into smaller subsets.When the subset limit after branching exceeds the known feasible solution value,the branching is terminated to reduce the search range until the feasible solution is found.

In the numerical expression in the following table, define “NumLPs”represents the number of traversal nodes during the iteration of a specific algorithm,“F”indicates whether the specific method is feasible for solving the example,including the following three cases.“Feasible”means that an integer solution can be obtained by calculating an instance,otherwise it is represented by“Infeasible”.We use the fixed-point iterative method and branch-and-bound method introduced in this study to solve the problems ofm= 3,n= 20,m= 4,n= 30 andm= 5,n= 40.Table 1 illustrates the numerical outcomes ofm=3 andn=20 problem,Table 2 presents the outcomes ofm=4 andn=30 problem,and the outcomes of the problem withm= 5 andn= 40 are shown in Table 3.Figs.4–6 show the comparison of the number of nodes traversed by the fixed-point iterative method and the branchand-bound method for the three dimensional problems mentioned above,respectively,which show the computational effects and performance of the two algorithms more clearly and intuitively.

The three groups of experimental values mentioned above show that the method studied in this paper can effectively solve the problem of DT.In the case ofm= 3,n= 20 in Table 1, there is no significant difference between the fixed-point iterative method and the branch and bound method in solving this problem,and the results can be obtained within 5 s,although the number of nodes traversed by the branch and bound method in the iteration is slightly larger than that of the fixed-point iterative method.Table 2 shows that, when the dimension of the problem increases tom= 4,n= 30, the fixed-point iterative method can generally end the solution process of the problem in 200 s and obtain the feasible solution.The number of nodes traversed by the branch-and-bound method in computing the feasible solution is about twice or more than that of the fixed-point iteration method because the dimension is so high at this point.In the two dimensional cases mentioned above, branch-andbound and fixed-point iterative methods are effective in computing feasible solutions to the problem without considering the number of traversed nodes.Finally,whenm=5,n=40,as shown in Table 3,the branch-and-bound method does not solve the problem in some cases.If the fixed-point iterative method is feasible without considering the length of time,the time to solve this problem may be longer or even several hours due to the different parameters selected in different examples.Therefore, we choose to terminate the iteration of the algorithm for solving this instance and mark it“Infeasible”.For the fixed-point iterative method,if the time consumption is ignored,it can solve the reconstruction model of DT problem with any dimension.

Table 1: m=3,n=20

Table 2: m=4,n=30

Table 2 (continued)Prob. Fixed-point Branch-and-bound F NumLPs F NumLPs 9 Feasible 220 Feasible 837 10 Feasible 340 Feasible 796

Table 3: m=5,n=40

From the data of the three charts mentioned above with increasing dimension,the number of nodes traversed in the process of obtaining the result through fixed-point iteration is about 40%higher than that of branch and bound method.With the continuous increase in dimension,the advantage of fixedpoint iteration method in traversing less nodes is more prominent,which also verifies the efficiency of this algorithm.

The comparison and analysis of the three tables mentioned above, indicate that solving the problem becomes increasingly difficult when the dimension of the problem increases due to the limited computing power of a single computer.If the fixed-point iterative method can be in a distributed implementation and multiple computers are used to solve this problem at the same time, then even higher dimensional problems can be solved with relatively less time consumption.Therefore,our future work direction will be to realize the distributed computing of fixed-point iterative method.

Figure 4:Comparison of the number of traversal nodes with m=3 and n=20

Figure 5:Comparison of the number of traversal nodes with m=4 and n=30

After the two methods mentioned above are compared,we change the fixed-point iterative method in simple steps to change the initial point of its iteration for forming an inverse iterative process.For the three types of problems in Tables 1–3,the fixed-point iterative method before and after changing the initial point is used to calculate each example.The number of iterations required in the process of calculating whether the problem has a solution is recorded in Fig.7,the data show that changing the initial point of iteration will not affect the feasibility and effectiveness of the fixed-point iterative method proposed in this study.

Figure 6:Comparison of the number of traversal nodes with m=5 and n=40

Figure 7:Iterative effect of different initial points

6 Conclusions and Future Prospects

This study presents a new method that can effectively solve the reconstruction problem of the DT model.The solution procedure is divided into two steps.First,the reconstruction problem is changed into a polyhedron judgment problem on the basis of LLL-based reduction algorithm.Second,the fixed point iteration method of integer programming is used to solve the polyhedron judgment problem in the previous step.The calculation results of specific examples are compared and analyzed between this method and branch and bound method.The results show that the fixed point iteration method has more advantages in solving this kind of problem which is embodied in the optimization of time and the number of traversal nodes.

However,the dimension of the problem examples involved in this study is not particularly high.Thus,the latent advantages of this method are not thoroughly reflected.With the continuous increase in the problem dimension, this method can be realized by distributed approach in the future.This solution process will reduce considerable time consumption through utilizing multiple slave computers working at the same time.

Funding Statement:This research was funded by the NSFC under Grant Nos.61803279, 71471091,62003231 and 51874205,in part by the Qing Lan Project of Jiangsu,in part by the China Postdoctoral Science Foundation under Grant Nos.2020M671596 and 2021M692369, in part by the Suzhou Science and Technology Development Plan Project (Key Industry Technology Innovation) under Grant No.SYG202114,in part by the Natural Science Foundation of Jiangsu Province under Grant No.BK20200989,and Postdoctoral Research Funding Program of Jiangsu Province.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.