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.

1. 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 [58]. In this paper, we propose a dynamical method based on the continuous-time optimization approach with nonlinear differential equations [913] 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.

2. Problem Description

Let be a set of positive real numbers, with projection and projection operator both given in advance. Define a set of labels , , that are the gray values [15]. Assume that Define also a vector and the corresponding matrix where is a identity matrix, indicates the transpose of a vector or matrix, and is the Kronecker product. A pixel vector ,  , is given by

With the above preliminaries, the discrete tomography described in this paper involves solving the following equation for unknown vector .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 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 . Then, the matrix elements are given as that is,wherefor 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 is a diagonal matrix whose diagonal elements are those of the vector . The matrix is written asNote that, for the true vector , the definition of guarantees thatTherefore, is definitely included as an equilibrium point of (9).

We can rewrite the dynamical system in (9) aswhere , and If we provide a matrix, then the following equation holds.

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 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.

3. Theoretical Analysis

We rewrite (9) asand 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 noise-free.

Proposition 1. If we choose initial value in the switched dynamical system in (9), then the solution stays in for all .

Proof. As the system can be written as , we see that, on the subspace where or , the solution satisfies 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 of the system in (9) with initial value being in for all .

The Jacobian or the derivative of with respect to isWe 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 . Besides the true equilibrium, other false equilibria exist in , described bywhile satisfying . Examples of arefor andfor . Next, we consider the sets for the true and false equilibria, respectively,

Proposition 3. If there exists an equilibrium in of (9), then it is locally half-stable.

Proof. When the equilibrium is in , which is a subset of , the Jacobian at point is given by because and . A diagonal matrix with nonpositive diagonal elements has nonpositive eigenvalues. This implies that the equilibrium is half-stable. However, for , we have the Jacobian at asThe 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 has the same eigenvalues as the matrix , which is a positive semidefinite matrix, where denotes a diagonal matrix satisfying . Then, all of the eigenvalues of the Jacobian are nonpositive and, therefore, is a half-stable equilibrium.

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

Some elements of are in , and the relationship between and isif and are nonempty sets.

Proposition 4. If contains an equilibrium of (9), then it is a saddle-type equilibrium.

Proof. The local stability of the equilibrium is determined by the eigenvalues of the Jacobian as From the definition of and , this can be rewritten asWe 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.

4. 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):whereBy multiplying or , 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.

5. Numerical Experiments

5.1. Simplest Numerical Experiment

We begin with the simplest possible example, that of a -pixel case. We set , , and , and defined asAccording to the above settings, we haveThe projection operator is given as

Now, we suppose a true solution asFigure 1(a) shows an example of a -pixel phantom image, and Figures 1(b) and 1(c) show true discrete images corresponding to and , respectively. Each pixel value , where , was determined by (3); that is, a pixel was expressed as a linear combination of unknowns and labels and . Given that we knew the true solution in advance, we could compute the projection data 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 and . These images are given asand 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 for any as .

We solved (9) by using the MATLAB function  ode113. Figure 2 shows the time responses of the orbits and for . 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 plane with two different initial values. The blue curve starting from 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 . First, the red curve runs along the separatrix and approaches ; however, it later turns to an equilibrium 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 . In other words, many possibilities for the system becoming trapped at an equilibrium in were excluded. Theoretically, the saddle separatrices split the 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 , we found that, for , the system converged to by . The resolution used for was . Even if the initial value takes the same value for each element of , a wide interval is available for the domain of attraction, which makes it reasonable for practical use.

5.2. (64 × 64)-Pixel Image Reconstruction

We prepared a -pixel digital phantom based on the Shepp-Logan model [16]; see Figure 4(a). Figure 4(b) illustrates the pixel-value 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 every . Thus, there were 90 directions in total; that is, . 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 asWe used the solver  ode113  to simulate (9). The initial value was .

From the above setup, we obtained the segmented images shown in Figures 5(a) and 5(b) according to the label values by computing and , respectively. This figure was computed by using various initial values . The pixel values in Figures 5(a) and 5(b) are either black or white; thus, the solutions belong to . Indeed, these solutions are in the true equilibrium solution set .

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 -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.

5.3. (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 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 -norm was 11. This relatively small value shows that self-adjusting labeling can be realized.

6. 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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.