# Multivalued Discrete Tomography Using Dynamical System That Describes Competition

Authors: Takeshi Kojima; Tetsushi Ueta; Tetsuya Yoshinaga
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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.


---

## Body

## 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 [5–8]. In this paper, we propose a dynamical method based on the continuous-time optimization approach with nonlinear differential equations [9–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.


## 2. Problem Description

LetR+ be a set of positive real numbers, with projection y∈R+I and projection operator A∈R+I×J both given in advance. Define a set L≔{g1,g2,…,gL} of labels gj∈(0,1], j=1,2,…,L, that are the gray values [15]. Assume that (1)0<g1<g2<⋯<gL.Define also a vector g≔(g1,g2,…,gL)⊤ and the corresponding matrix (2)G≔UJ⊗g⊤∈0,1J×JL,where UJ is a J×J identity matrix, ⊤ indicates the transpose of a vector or matrix, and ⊗ is the Kronecker product. A pixel vector c≔(c1,c2,…,cJL)⊤,  cj∈[0,1], is given by(3)c=Gx.With the above preliminaries, the discrete tomography described in this paper involves solving the following equation for unknown vectorx∈[0,1]JL.(4)y=AGx.Note that Gx represents the gray values in a reconstructed image. Ideally, each element of x should be a binary number, but, here, we assume it is a real number in the interval [0,1] to accommodate cases in which gj is given incorrectly or the inverse problem is well posed.If the problem is consistent, a true solution of (4) is denoted as e∈{0,1}JL. Then, the matrix elements are given as (5)djl≔ej-1L+l,j=1,2,…,J,l=1,2,…,L;that is,(6)D≔d11d12⋯d1Ld21d22⋯d2L⋮⋮⋱⋮dJ1dJ2⋯dJL∈RJ×L,where(7)∑l=1Ldjl=1for any jth row. However, if the corresponding pixel is in the background of an image, we have(8)∑l=1Ldjl=0.To solve (4), we utilize a dynamical system approach; that is, we rewrite the problem as an initial value problem of a differential equation:(9)dxdt=XUJL-XAG⊤y-AGx-Ψx,x0=x0,where X≔diag(x) is a diagonal matrix whose diagonal elements are those of the vector x. The matrix Ψ is written as(10)Ψ≔1LUJ⊗gg⊤-diag2⁡g.Note that, for the true vector e, the definition of Ψ guarantees that(11)diag⁡eΨe=0.Therefore, e is definitely included as an equilibrium point of (9).We can rewrite the dynamical system in (9) as(12)dwjldt=wjl1-wjlAG⊤j-1L+ly-AGx-glL∑k=1k≠lLgkwjk,where j=1,2,…,J,l=1,2,…,L, and (13)wjl≔xj-1L+l.If we provide a matrix, (14)W≔w11w12⋯w1Lw21w22⋯w2L⋮⋮⋱⋮wJ1wJ2⋯wJL∈RJ×L,then the following equation holds. (15)Gx=Wg.In (9), without Ψx, 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 X(U-X) is appended in anticipation of the solution x(t) 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 (-Ψx), is inspired by the generalized Lotka-Volterra equation [14] to ensure (16)∑l=1Lwjlt⟶1,t⟶∞;namely, for some l,(17)∑k=1k≠lLwjkt⟶0,wjlt⟶1,t⟶∞in (12) for j such that the condition in (7) is satisfied.


## 3. Theoretical Analysis

We rewrite (9) as(18)dxdt=fxand assume that A and y 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 valuex0∈(0,1)JL in the switched dynamical system in (9), then the solution ϕ(t,x0) stays in (0,1)JL for all t∈R+.Proof.
As the system can be written asdxj/dt=fj(x), we see that, on the subspace where xj=0 or xj=1, the solution satisfies dϕj/dt≡0 for any j. 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 ϕ(t,x0) of the system in (9) with initial value x0∈(0,1)JL being in (0,1)JL for all t∈R+.The Jacobian or the derivative off with respect to x is(19)∂f∂xx=-XUJL-XAG⊤AG+Ψ+UJL-2Xdiag⁡AG⊤y-AGx-Ψx.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 u, are, respectively, (20)∂f∂x0=diag⁡AG⊤y=diag⁡AG⊤AGe,∂f∂xu=diag⁡AG⊤-y+AGu+Ψu=diag⁡AG⊤AGu-e+Ψu.We see that all of the eigenvalues of each matrix are nonnegative, and accordingly, both equilibria are unstable.Let us define the set(21)S≔s∈0,1JL:y-AGs=0,diag⁡sΨs=0.Note that the exact equilibrium e belongs to S. Besides the true equilibrium, other false equilibria exist in S, described by(22)UJ⊗He,while satisfying g⊤H=g⊤. Examples of H are(23)00g1g21for L=2 and(24)000000g1g3g2g31,000g1g2000g2g31,000g1g210001for L=3. Next, we consider the sets for the true and false equilibria, respectively,(25)T≔S∩0,1JL,F≔S∖T.Proposition 3.
If there exists an equilibrium inS of (9), then it is locally half-stable.Proof.
When the equilibriume is in T, which is a subset of S, the Jacobian at point e is given by (26)∂f∂xe=-UJL-2diag⁡ediag⁡Ψe=-diag⁡Ψebecause diag⁡e(UJL-diag⁡e)=0 and diag⁡eΨe=0. A diagonal matrix with nonpositive diagonal elements has nonpositive eigenvalues. This implies that the equilibrium e∈T is half-stable. However, for s∈F, we have the Jacobian at s as(27)∂f∂xs=-diag⁡sUJL-diag⁡sAG⊤AG-diag⁡sUJL-diag⁡sΨ-diag⁡Ψs.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⁡s(UJL-diag⁡s)(AG)⊤AG has the same eigenvalues as the matrix S(AG)⊤(AG)S, which is a positive semidefinite matrix, where S denotes a diagonal matrix satisfying S2=diag⁡s(UJL-diag⁡s). Then, all of the eigenvalues of the Jacobian are nonpositive and, therefore, s 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(28)Z≔z∈0,1JL∖S:AG⊤y-AG⊤AG+Ψz=0,Z¯≔z¯∈0,1JL∖S∪0JL∪1JL:diag⁡z¯UJL-diag⁡z¯AG⊤y-AG⊤AG+Ψz¯=0.Some elements ofz¯∈Z¯ are in {0,1}, and the relationship between z∈Z and z¯∈Z¯ is(29)diag⁡z¯UJL-diag⁡z¯AG⊤AG+Ψz-z¯=0if Z and Z¯ are nonempty sets.Proposition 4.
IfZ contains an equilibrium of (9), then it is a saddle-type equilibrium.Proof.
The local stability of the equilibriumz∈Z is determined by the eigenvalues of the Jacobian as (30)∂f∂xz=-diag⁡zUJL-diag⁡zAG⊤AG+Ψ.From the definition of G and Ψ, this can be rewritten as(31)∂f∂xz=-diag⁡zUJL-diag⁡zA⊤A+1LUJ⊗gg⊤+diag⁡zUJL-diag⁡z1LUJ⊗diag2⁡g.We see that the matrix of the first term has rank J 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 ∂f/∂x(z) or equivalently the trace of the matrix -diag(z)(UJL-diag(z))(AG)⊤AG, 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 y=AGx 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):(32)dxdt=XUJL-XαtAG⊤y-AGx-1-αtΨx,where(33)αt=exp⁡-tτ.By multiplying (α(t)) or (1-α(t)), the effects of the term (AG)⊤(y-AGx) or the term (-Ψx) are emphasized or restrained by the parameter τ.In the system given by (32), at early times t, the orbit is affected by the term (AG)⊤(y-AGx); thus, Wg approaches the nondiscrete reconstructed image Dg. This effect is gradually restrained as the effect of the term (-Ψx) becomes dominant as t 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(2×2)-pixel case. We set J=4, I=6, and L=2, and defined g as(34)g=0.51.According to the above settings, we have(35)G=0.51000000000.51000000000.51000000000.51,Ψ=00.250000000.2500000000000.250000000.2500000000000.250000000.2500000000000.250000000.250.The projection operator is given as(36)A=101001011001001111000110.Now, we suppose a true solution as(37)e1e2e3e4e5e6e7e8=10010100.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 g1 and g2, respectively. Each pixel value ci, where i=1,2,3,4, was determined by (3); that is, a pixel was expressed as a linear combination of unknowns x=(x1,x2,…,x8)⊤ and labels g1 and g2. Given that we knew the true solution e=(e1,e2,…,e8)⊤ in advance, we could compute the projection data y=(y1,y2,…,y6)⊤=(1.5,1,0.5,1,1.5,2)⊤ by evaluating y=AGe. Let us describe the problem in this paper again, with A and G given. We solve for unknowns x from the projection y given by a measurement. Consequently, we obtain discrete reconstructed images corresponding to labels g1 and g2. These images are given as(38)x1x3x5x7,x2x4x6x8,and the image composed of them is expressed as(39)C=c1c2c3c4=g1x1+g2x2g1x3+g2x4g1x5+g2x6g1x7+g2x8.Figure 1
(
2
×
2
)-pixel phantom. White, gray, and black pixels are unity, 0.5, and zero values, respectively.
(a)
(b)
(c)The trajectoryx 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 x=e. 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 xj0=u for any j as x∗0=u.We solved (9) by using the MATLAB function  ode113. Figure 2 shows the time responses of the orbits x1(t) and x2(t) for x∗0=0.85. The orbit transiently approached and left the saddle z and asymptotically converged toward the stable equilibrium e.Figure 2
Time response ofx1 and x2 starting from x∗0=0.85 in subdivisions according to pixels.Let us discuss the initial value dependency of (9) for this example. Figure 3 shows a phase portrait in the x1-x2 plane with two different initial values. The blue curve starting from x∗0=0.77 approaches the saddle equilibrium point z once before finally converging to the true solution e. The green and magenta dotted lines are the stable and unstable manifolds of z, respectively, and the red curve shows an orbit starting from x∗0=0.72. First, the red curve runs along the separatrix and approaches z; however, it later turns to an equilibrium s∈F corresponding to a false solution. The location of s is(40)s1s2s3s4s5s6s7s8=00.5010100.Figure 3
Equilibria and attractors projected onx1-x2 phase plane.From Proposition4, the term (-Ψx) 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)JL 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 e 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 x∗0=x∗, we found that, for 0.7443<x∗<0.9950, the system converged to e by t=400. The resolution used for x∗ was 10-4. Even if the initial value takes the same value for each element of x0, 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(64×64)-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 180∘ every 2∘. Thus, there were 90 directions in total; that is, I=95×90. The projection operator A 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(41)g=0.51.We used the solver  ode113  to simulate (9). The initial value was x∗0=0.85.Figure 4
(a) Phantom image, (b) pixel value of each segment, (c) sinogram, and (d) reconstructed image obtained by FBP.
(a)
(b)
(c)
(d)From the above setup, we obtained the segmented images shown in Figures5(a) and 5(b) according to the label values by computing W(t)(1,0)⊤ and W(t)(0,1)⊤, respectively. This figure was computed by using various initial values 0<x∗0<1. The pixel values in Figures 5(a) and 5(b) are either black or white; thus, the solutions x belong to {0,1}JL. Indeed, these solutions are in the true equilibrium solution set T.Figure 5
Reconstructed image. (a) Image ofW(t)(1,0)⊤, (b) image of W(t)(0,1)⊤, and (c) image of W(t)g with x∗0=0.85 at t=10,000.
(a)
(b)
(c)Figure5(c) shows a composited image obtained by computing W(t)g. 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 L1-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.Figure 6
Density profiles. Gray values of (a) phantom image and (b) images with proposed method and FBP. Green, red, and blue lines indicate gray values of images shown in Figures4(a), 4(d), and 5(c), respectively.
(a)
(b)


### 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 Figure7(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 Dg image corresponding to the given label set g=(0.5,1)⊤ is shown in Figure 7(e). This figure is the objective image for the composited images generated by the discrete tomography.Figure 7
(a) Phantom image, (b) pixel value of each segment, (c) sinogram, (d) reconstructed image obtained by FBP, and (e) objective composited image.
(a)
(b)
(c)
(d)
(e)The results are shown in Figures8(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 L1-norm was 11. This relatively small value shows that self-adjusting labeling can be realized.Figure 8
Reconstructed image. (a) Image ofW(t)(1,0)⊤, (b) image of W(t)(0,1)⊤, and (c) image of W(t)g with x∗0=0.85 and τ=1,000 at t=400,000.
(a)
(b)
(c)


## 5.1. Simplest Numerical Experiment

We begin with the simplest possible example, that of a(2×2)-pixel case. We set J=4, I=6, and L=2, and defined g as(34)g=0.51.According to the above settings, we have(35)G=0.51000000000.51000000000.51000000000.51,Ψ=00.250000000.2500000000000.250000000.2500000000000.250000000.2500000000000.250000000.250.The projection operator is given as(36)A=101001011001001111000110.Now, we suppose a true solution as(37)e1e2e3e4e5e6e7e8=10010100.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 g1 and g2, respectively. Each pixel value ci, where i=1,2,3,4, was determined by (3); that is, a pixel was expressed as a linear combination of unknowns x=(x1,x2,…,x8)⊤ and labels g1 and g2. Given that we knew the true solution e=(e1,e2,…,e8)⊤ in advance, we could compute the projection data y=(y1,y2,…,y6)⊤=(1.5,1,0.5,1,1.5,2)⊤ by evaluating y=AGe. Let us describe the problem in this paper again, with A and G given. We solve for unknowns x from the projection y given by a measurement. Consequently, we obtain discrete reconstructed images corresponding to labels g1 and g2. These images are given as(38)x1x3x5x7,x2x4x6x8,and the image composed of them is expressed as(39)C=c1c2c3c4=g1x1+g2x2g1x3+g2x4g1x5+g2x6g1x7+g2x8.Figure 1
(
2
×
2
)-pixel phantom. White, gray, and black pixels are unity, 0.5, and zero values, respectively.
(a)
(b)
(c)The trajectoryx 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 x=e. 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 xj0=u for any j as x∗0=u.We solved (9) by using the MATLAB function  ode113. Figure 2 shows the time responses of the orbits x1(t) and x2(t) for x∗0=0.85. The orbit transiently approached and left the saddle z and asymptotically converged toward the stable equilibrium e.Figure 2
Time response ofx1 and x2 starting from x∗0=0.85 in subdivisions according to pixels.Let us discuss the initial value dependency of (9) for this example. Figure 3 shows a phase portrait in the x1-x2 plane with two different initial values. The blue curve starting from x∗0=0.77 approaches the saddle equilibrium point z once before finally converging to the true solution e. The green and magenta dotted lines are the stable and unstable manifolds of z, respectively, and the red curve shows an orbit starting from x∗0=0.72. First, the red curve runs along the separatrix and approaches z; however, it later turns to an equilibrium s∈F corresponding to a false solution. The location of s is(40)s1s2s3s4s5s6s7s8=00.5010100.Figure 3
Equilibria and attractors projected onx1-x2 phase plane.From Proposition4, the term (-Ψx) 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)JL 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 e 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 x∗0=x∗, we found that, for 0.7443<x∗<0.9950, the system converged to e by t=400. The resolution used for x∗ was 10-4. Even if the initial value takes the same value for each element of x0, 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(64×64)-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 180∘ every 2∘. Thus, there were 90 directions in total; that is, I=95×90. The projection operator A 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(41)g=0.51.We used the solver  ode113  to simulate (9). The initial value was x∗0=0.85.Figure 4
(a) Phantom image, (b) pixel value of each segment, (c) sinogram, and (d) reconstructed image obtained by FBP.
(a)
(b)
(c)
(d)From the above setup, we obtained the segmented images shown in Figures5(a) and 5(b) according to the label values by computing W(t)(1,0)⊤ and W(t)(0,1)⊤, respectively. This figure was computed by using various initial values 0<x∗0<1. The pixel values in Figures 5(a) and 5(b) are either black or white; thus, the solutions x belong to {0,1}JL. Indeed, these solutions are in the true equilibrium solution set T.Figure 5
Reconstructed image. (a) Image ofW(t)(1,0)⊤, (b) image of W(t)(0,1)⊤, and (c) image of W(t)g with x∗0=0.85 at t=10,000.
(a)
(b)
(c)Figure5(c) shows a composited image obtained by computing W(t)g. 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 L1-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.Figure 6
Density profiles. Gray values of (a) phantom image and (b) images with proposed method and FBP. Green, red, and blue lines indicate gray values of images shown in Figures4(a), 4(d), and 5(c), respectively.
(a)
(b)


## 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 Figure7(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 Dg image corresponding to the given label set g=(0.5,1)⊤ is shown in Figure 7(e). This figure is the objective image for the composited images generated by the discrete tomography.Figure 7
(a) Phantom image, (b) pixel value of each segment, (c) sinogram, (d) reconstructed image obtained by FBP, and (e) objective composited image.
(a)
(b)
(c)
(d)
(e)The results are shown in Figures8(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 L1-norm was 11. This relatively small value shows that self-adjusting labeling can be realized.Figure 8
Reconstructed image. (a) Image ofW(t)(1,0)⊤, (b) image of W(t)(0,1)⊤, and (c) image of W(t)g with x∗0=0.85 and τ=1,000 at t=400,000.
(a)
(b)
(c)


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



---

*Source: 8160354-2017-11-12.xml*