Multivalued Discrete Tomography Using Dynamical System That Describes Competition

Multivalued discrete tomography involves reconstructing images composed of three or more gray levels from projections. We propose a method based on the continuous-time optimization approach with a nonlinear dynamical system that effectively utilizes competition dynamics to solve the problem of multivalued discrete tomography. We perform theoretical analysis to understand how the system obtains the desired multivalued reconstructed image. Numerical experiments illustrate that the proposed method also works well when the number of pixels is comparatively high even if the exact labels are unknown.


Introduction
Multivalued (or nonbinary) discrete tomography involves the reconstruction of images composed of three or more gray levels from projections.Compared with computed tomography, it is possible to reduce the number of projections by using prior knowledge about a set of gray levels.This is important for medical use as it is the basis for identifying characteristic regions in tomographic images [1,2].Conventional methods for discrete tomography include the iterative reconstruction method involving an iterative algorithm and image segmentation [3], an optimization algorithm based on minimizing an energy function to discretize multiple intensity values [4], and various other methods [5][6][7][8].In this paper, we propose a dynamical method based on the continuous-time optimization approach with nonlinear differential equations [9][10][11][12][13] that are capable of obtaining a desired tomographic image through convergence to a limit set of the differential equations.Our method utilizes the competition dynamics of generalized Lotka-Volterra systems [14] to solve the problem of multivalued discrete tomography.A nonlinear term that conducts the competitive behavior of a solution is incorporated into differential equations to ensure that the solution orbit starting from an appropriate initial value converges satisfactorily to the desired solution.
We propose two differential equations to represent an autonomous system and a nonautonomous system that have similar structures.For the autonomous system, it has been proven theoretically that the stable equilibria corresponding to the ideal solution and the undesired solution coexist and that a saddle-type equilibrium exists that plays an important role in the behavior of the solutions.We investigate the mechanism behind this behavior through a numerical example.The results of numerical experiments show that the proposed method works well even if the number of pixels is comparatively high.Further numerical experiments show that the nonautonomous system can be applied in cases in which the exact label set is not given.

𝑦 = 𝐴𝐺𝑥.
(4) Note that  represents the gray values in a reconstructed image.Ideally, each element of  should be a binary number, but, here, we assume it is a real number in the interval [0, 1] to accommodate cases in which   is given incorrectly or the inverse problem is well posed.
If the problem is consistent, a true solution of ( 4) is denoted as  ∈ {0, 1}  .Then, the matrix elements are given as where for any th row.However, if the corresponding pixel is in the background of an image, we have To solve (4), we utilize a dynamical system approach; that is, we rewrite the problem as an initial value problem of a differential equation: where  fl diag() is a diagonal matrix whose diagonal elements are those of the vector .The matrix Ψ is written as Note that, for the true vector , the definition of Ψ guarantees that diag () Ψ = 0.
In (9), without Ψ, the dynamics is based on gradient systems proposed for binary tomography inverse problems in [9,11].In the former reference, a dynamical system is provided whose vector field resembles that of a logistic equation, and the convergence of solutions is demonstrated theoretically.In the latter reference, a further term ( − ) is appended in anticipation of the solution () wandering inside of the subspace (0, 1) and converging either to the true value zero or to unity.
In this paper, we treat multivalued discrete tomography as an extension of the binary tomography problems addressed in [9,11].Equation (9) shows that the proposed system, including the term (−Ψ), is inspired by the generalized Lotka-Volterra equation [14] to ensure namely, for some ℓ, in (12) for  such that the condition in ( 7) is satisfied.

Theoretical Analysis
We rewrite (9) as and assume that  and  are nonnegative.Equation ( 18) has equilibria that include the zero vector, the vector whose elements are all unity, and a constant nonzero vector that corresponds to the desired reconstructed image, assuming that the projection data are complete, consistent, and noisefree.
Proof.As the system can be written as   / =   (), we see that, on the subspace where   = 0 or   = 1, the solution satisfies   / ≡ 0 for any .Therefore, the subspace is invariant, and trajectories cannot pass through every invariant subspace, according to the uniqueness of solutions for the initial value problem.This leads to any solution (,  0 ) of the system in ( 9) with initial value  0 ∈ (0, 1)  being in (0, 1)  for all  ∈  + .
The Jacobian or the derivative of  with respect to  is We can prove propositions concerning local stability as follows.
Proposition 2. Each of the all-zeros and all-ones equilibria of ( 9) is locally unstable.
Proof.From (19), the Jacobian matrices at the all-zeros equilibrium and all-ones equilibrium, say , are, respectively, We see that all of the eigenvalues of each matrix are nonnegative, and accordingly, both equilibria are unstable.
Let us define the set Note that the exact equilibrium  belongs to S.Besides the true equilibrium, other false equilibria exist in S, described by while satisfying  ⊤  =  ⊤ .Examples of  are for  = 2 and ) , for  = 3. Next, we consider the sets for the true and false equilibria, respectively, Proposition 3. If there exists an equilibrium in S of ( 9), then it is locally half-stable.
Proof.When the equilibrium  is in T, which is a subset of S, the Jacobian at point  is given by because diag()(  − diag()) = 0 and diag()Ψ = 0.
A diagonal matrix with nonpositive diagonal elements has nonpositive eigenvalues.This implies that the equilibrium  ∈ T is half-stable.However, for  ∈ F, we have the Jacobian at  as The eigenvalues of the Jacobian are the sum of the eigenvalues of each of the three terms in (27); note that the second term has all-zero eigenvalues.The first term in (27) is a negative semidefinite matrix because diag()(  −diag())(() ⊤ ) has the same eigenvalues as the matrix () ⊤ (), which is a positive semidefinite matrix, where  denotes a diagonal matrix satisfying  2 = diag()(  − diag()).Then, all of the eigenvalues of the Jacobian are nonpositive and, therefore,  is a half-stable equilibrium.

Mathematical Problems in Engineering
Numerous saddle-type equilibria exist in the system, and these play an important role in separating trajectories that converge to true and false stable equilibria.We consider the two equilibrium sets (28) Some elements of  ∈ Z are in {0, 1}, and the relationship between  ∈ Z and if Z and Z are nonempty sets.Proposition 4. If Z contains an equilibrium of ( 9), then it is a saddle-type equilibrium.
Proof.The local stability of the equilibrium  ∈ Z is determined by the eigenvalues of the Jacobian as From the definition of  and Ψ, this can be rewritten as We see that the matrix of the first term has rank  and all its eigenvalues are nonpositive.The second term is a diagonal matrix of full rank with positive eigenvalues, so the Jacobian has positive eigenvalues.However, from (7), with Ψ having zero diagonal elements, the sum of the eigenvalues is the trace of /() or equivalently the trace of the matrix , which is negative.Therefore, the eigenvalues include both the positive and negative values, meaning that the equilibrium is a saddle.

Promotion of Distinction
Our proposed system in ( 9) can obtain a solution that resolves the tomographic inverse problem of  =  and satisfies the conditions in (7) or (8), when assumed that the exact label set is given.To relax the assumption and satisfy (7) or ( 8) even if no exact label set is given, that is, realize image segmentation based on the labels with a small range of gray values rounded to the nearest gray label, we propose a nonautonomous system that is an improvement of the system given by ( 9): where By multiplying (()) or (1 − ()), the effects of the term () ⊤ ( − ) or the term (−Ψ) are emphasized or restrained by the parameter .
In the system given by (32), at early times , the orbit is affected by the term () ⊤ ( − ); thus,  approaches the nondiscrete reconstructed image .This effect is gradually restrained as the effect of the term (−Ψ) becomes dominant as  grows.Consequently, the state variables from which a pixel value is structured begin to compete with each other.Therefore, one of the state variables is enforced to be nonzero, and the others are zeros.We refer to this effect as self-adjusting labeling.
According to the above settings, we have  = ( 0.5 1 0 0 0 0 0 0 0 0 0.5 1 0 0 0 0 0 0 0 0 0.5 1 0 0 0 0 0 0 0 0 0.5  The projection operator is given as Now, we suppose a true solution as Figure 1(a) shows an example of a (2 × 2)-pixel phantom image, and Figures 1(b) and 1(c) show true discrete images corresponding to  1 and  2 , respectively.Each pixel value   , where  = 1, 2, 3, 4, was determined by (3); that is, a pixel was expressed as a linear combination of unknowns  = ( 1 ,  2 , . . .,  8 ) ⊤ and labels  1 and  2 .Given that we knew the true solution  = ( 1 ,  2 , . . .,  8 ) ⊤ in advance, we could compute the projection data  = ( 1 ,  2 , . . .,  6 ) ⊤ = (1.5, 1, 0.5, 1, 1.5, 2) ⊤ by evaluating  = .Let us describe the problem in this paper again, with  and  given.We solve for unknowns  from the projection  given by a measurement.Consequently, we obtain discrete reconstructed images corresponding to labels  1 and  2 .These images are given as and the image composed of them is expressed as The trajectory  obeying (9) starts from arbitrary initial values, and if the orbit converges to a stable equilibrium point, we expect that the system will offer the true solution  = .Because no clues as to the initial values can be obtained a priori, all elements of the initial value are simply aligned with the same value.Let us denote  0  =  for any  as  0 * = .We solved (9) by using the MATLAB function ode113.Figure 2 shows the time responses of the orbits  1 () and  2 () for  0 * = 0.85.The orbit transiently approached and left the saddle  and asymptotically converged toward the stable equilibrium .
Let us discuss the initial value dependency of (9) for this example.Figure 3 shows a phase portrait in the  1 - 2 plane Mathematical Problems in Engineering with two different initial values.The blue curve starting from  0 * = 0.77 approaches the saddle equilibrium point  once before finally converging to the true solution .The green and magenta dotted lines are the stable and unstable manifolds of , respectively, and the red curve shows an orbit starting from  0 * = 0.72.First, the red curve runs along the separatrix and approaches ; however, it later turns to an equilibrium  ∈ F corresponding to a false solution.The location of  is From Proposition 4, the term (−Ψ) in ( 9) means that quite a few equilibria were turned into saddle-type ones that belonged in Z.In other words, many possibilities for the system becoming trapped at an equilibrium in Z were excluded.Theoretically, the saddle separatrices split the (0, 1)  space in (9) definitively into several domains of attraction.However, we could not find an appropriate set of initial values in the domain of attraction for  because the system was high dimensional.Instead, for this experiment, we checked the domain of attraction by brute force.By rewriting the initial value as  0 * =  * , we found that, for 0.7443 <  * < 0.9950, the system converged to  by  = 400.The resolution used for  * was 10 −4 .Even if the initial value takes the same value for each element of  0 , a wide interval is available for the domain of attraction, which makes it reasonable for practical use.

(64 × 64)-Pixel Image
Reconstruction.We prepared a (64 × 64)-pixel digital phantom based on the Shepp-Logan model [16]; see Figure 4(a).Figure 4(b) illustrates the pixelvalue distribution of each segment.We simulated a scanner that is equipped with 95 X-ray detectors per row and acquired parallel-beam projection data over 180 ∘ every 2 ∘ .Thus, there were 90 directions in total; that is,  = 95 × 90.The projection operator  was obtained by computing the discrete Radon transform.Figure 4(c) shows the corresponding sinogram.The reconstructed image obtained by the filtered back-projection (FBP) [17] with the Ramachandran-Lakshminarayanan filter from the sinogram is shown in Figure 4(d).Since the total number of projections was not sufficient, the reconstructed image did not match with the phantom image completely.The labels are defined as  = ( 0.5

1
) . (41) We used the solver ode113 to simulate (9).The initial value was  0 * = 0.85.From the above setup, we obtained the segmented images shown in Figures 5(a   It was clarified that the edges of our proposed method were sharper than the FBP's.In fact, the total difference with 1-norm between Figures 4(a) and 4(d) was 87; in comparison, the difference between Figures 4(a) and 5(c) was 1.3.Therefore, our proposed method generated a composited image, providing evidence that discrete tomography can produce accurate results.

(64 × 64)-Pixel Image Reconstruction with Self-Adjusting
Labeling.In the previous experiments, we assumed that the exact label set was given.The proposed nonautonomous system in (32), which is without this assumption, is aimed at reconstructing segmented binary images on the basis of given labels.Namely, we expect the system to automatically round a pixel that is not listed in the label set distinguished to the nearest label.Instead of using the system proposed in (9), we employed the system proposed in (9) to show the reconstruction result, wherein some pixel values do not match any element in the label set.
Let us provide a phantom that contains four different pixel values: 0.5, 0.6, 0.9, and 1; see Figure 7(a).Figure 7(b) illustrates the distribution of the pixel value of each segment.The parameter  was 1,000.The other conditions remained unchanged from those in Section 5.2. Figure 7(c) shows the corresponding sinogram.The reconstructed image obtained by the FBP from the sinogram is shown in Figure 7(d).The  image corresponding to the given label set  = (0.5, 1) ⊤ is shown in Figure 7(e).This figure is the objective image for the composited images generated by the discrete tomography.
The results are shown in Figures 8(a) and 8(b), where it was confirmed that the segmented images were obtained as intended by the nonautonomous system (32).As expected, in Figure 8(c), the gray values 0.6 and 0.9 were rounded to labels 0.5 and 1, respectively.When we compared the composited image in Figure 8(c) with the objective image composited in Figure 7(e), the difference with 1-norm was 11.This relatively small value shows that self-adjusting labeling can be realized.

Conclusion
We proposed a novel method for solving the problem of multivalued discrete tomography.Our method is based on the initial value problem of a nonlinear differential equation, which is inspired by the Lotka-Volterra competitive activity that enforces exclusivity among the state variables from which pixel values are constructed.
We proved the stability of all equilibria when the tomographic inverse problem was well posed.The equilibrium corresponding to the desired reconstructed image was stable; however, other false stable equilibria corresponding to undesired images coexisted.Therefore, a solution orbit that converged to a true or false equilibrium was determined by the initial value.
From the numerical experiments, we observed that a solution starting from the same uniform initial value converged to the true equilibrium, regardless of the pattern or the size of an image.Moreover, we proposed a modified system that is aimed at realizing self-adjusting labeling by adding a nonautonomous term.We confirmed that the nonautonomous system automatically classifies pixels that are not listed in the label set distinguished to the nearest label.

1 .
) and 5(b) according to the label values by computing ()(1, 0) ⊤ and ()(0, 1) ⊤ , respectively.This figure was computed by using various initial values 0 <  0 * < The pixel values in Figures5(a) and 5(b) are either black or white; thus, the solutions  belong to {0, 1}  .Indeed, these solutions are in the true equilibrium solution set T.

Figure 5 (
Figure 5(c) shows a composited image obtained by computing ().For comparison, let us show the density profile of the 26th row in each image.Figure 6(a) shows the density profile in the phantom image Figure 4(a), while, in Figure 6(b), the red and blue lines show the density profiles of the corresponding rows in Figures 4(d) and 5(c), respectively.It was clarified that the edges of our proposed method were sharper than the FBP's.In fact, the total difference with 1-norm between Figures4(a) and 4(d) was 87; in comparison, the difference between Figures4(a) and 5(c) was 1.3.Therefore, our proposed method generated a composited image, providing evidence that discrete tomography can produce accurate results.