Comparison of Sinogram-based Iterative Reconstruction with Compressed Sensing Techniques in X-ray CT

First, we review the sinogram-based iterative reconstruction (SbIR) algorithm [21]. Let $\mu_{1}$ be the density matrix (of size $n_{x}$ lines by $n_{y}$ columns) to be reconstructed from sinogram $y$, $n_{d}$ the number of detectors and $n_{v}$ the number of views. If $M$ is $n_{d}*n_{v}$ and $N$ is $n_{x}*n_{y}$, then the system matrix $A$ is of size $M$ lines by $N$ columns and the sinogram $y$ is a vector of size $M$. If we write $\mu_{1}$ in vector format then $\mu_{1}$ is a vector of size $N$; the reconstruction problem is to find $\mu_{1}$ such that

$$A*\mu_{1}=y.$$

There are two steps in the reconstruction process of the SbIR algorithm:

initialization:

first we initialize the reconstruction $\mu_{1}$ with an approximation derived from the sinogram $y$ as follows:

1. 1.1

   $\mu_{1}$ $\leftarrow$ $A^{\prime}$ $*$ ($y$ $./$ $\alpha$)

2. 1.2

   $\mu_{1}$ $\leftarrow$ $\mu_{1}$ $./$ $\beta$

iterations:

at each iteration, we update the reconstruction $\mu_{1}$ based on the difference between the sinogram $y$ and the current situation in $\mu_{1}$:

1. 2.1

   $\tilde{y}$ $\leftarrow$ $A$ $*$ $\mu_{1}$

2. 2.2

   $\mu_{2}$ $\leftarrow$ diag($\mu_{1}$) $*$ $A^{\prime}$ $*$ ($y$ $./$ $\tilde{y}$)

3. 2.3

   $\mu_{1}$ $\leftarrow$ $\mu_{2}$ $./$ $\beta$

where:

• •

  $A^{\prime}$ is the transpose of $A$,

• •

  $./$ is the element-by-element division,

• •

  $\tilde{y}$ is the situation in the reconstruction $\mu_{1}$ at the current iteration, and is a vector of size $M$ as the sinogram $y$,

• •

  $\alpha=\sum_{j}A[i,j]$ is of size $M$,

• •

  $\beta=\sum_{i}A[i,j]$ is of size $N$;

• •

  each element of $\mu_{1}$ corresponds to a square of size 1.0 by 1.0 in the experiment;

• •

  we also consider that $y$ has been obtained as area integrals, at a Detector,View pair meaning that each of the elements that are traversed is traversed for an area $\leq$ 1.0; we could consider line integrals as well, but using area integrals seem to help in the understanding of the algorithm.

The algorithm is very natural and its workings can be explained as follows; first let’s explain the initialization (when looking at the steps 1.1, 1.2 above, look only at the first element of $\mu_{1}$ which is $\mu_{1}[1]$):

• •

  in step 1.1, each of the $M$ elements of $y$ $./$ $\alpha$ is the average sinogram measurement per unit of area traversed by the X-ray beam at the respective Detector,View pair; thus each element of $y$ $./$ $\alpha$ can be viewed as an approximation of each of the elements of $\mu_{1}$ that are traversed by the X-ray beam at the respective Detector,View pair;

• •

  so, in step 1.1 the first element of $A^{\prime}$ $*$ ($y$ $./$ $\alpha$) is the sum of all approximations of $\mu_{1}[1]$ (from all Detector,View pairs in which $\mu_{1}[1]$ is traversed);

• •

  in step 1.2 $\mu_{1}[1]$ is divided by $\beta[1]$, where $\beta[1]$ is the sum of all non-zero elements in the first line of $A^{\prime}$; this is natural, because each of the approximations of $\mu_{1}[1]$ that come from step 1.1 are for unit of area 1.0, but for some Detector,View pairs in which $\mu_{1}[1]$ is traversed it might be traversed only for some portion of its total area 1.0.

Given this explanation of the initialization step, a better writing of the algorithm would be:

initialization:

first we initialize the reconstruction $\mu_{1}$ with an approximation derived from the sinogram $y$ as follows:

1. 1.1

   $\mu_{1}$ $\leftarrow$ $B$ $*$ ($y$ $./$ $\alpha$)

iterations:

at each iteration, we update the reconstruction $\mu_{1}$ based on the difference between the sinogram $y$ and the current situation in $\mu_{1}$:

1. 2.1

   $\tilde{y}$ $\leftarrow$ $A$ $*$ $\mu_{1}$

2. 2.2

   $\mu_{1}$ $\leftarrow$ diag($\mu_{1}$) $*$ $B$ $*$ ($y$ $./$ $\tilde{y}$)

where:

• •

  $B$ is a matrix of size $N$ lines by $M$ columns, whose $i$-th line is obtained from the $i$-th line of $A^{\prime}$ by dividing each element of the $i$-th line of $A^{\prime}$ to $\beta[i]$.

Given this re-writing of the algorithm, the workings can be explained as follows (again look only at $\mu_{1}[1]$):

• •

  in step 1.1, $\mu_{1}[1]$ is obtained (as a weighted sum) by multiplying the first line of $B$ with the vector of approximations $y$ $./$ $\alpha$; precisely, consider that there are $k$ ($k\leq{M}$) non-zero elements in the first line of $B$, let these be $B[1,j_{1}]$, $B[1,j_{2}]$,…, $B[1,j_{k}]$ where $1\leq{j_{1}}\leq{j_{2}}\leq...\leq{j_{k}}\leq{M}$; each of $B[1,j_{1}]$, $B[1,j_{2}]$,…, $B[1,j_{k}]$ corresponds to one of those Detector,View pairs that traverse $\mu_{1}[1]$, and represents the fraction (weight) from the total sum traversed $\beta[1]$; so we have

  $$B[1,j_{1}]+B[1,j_{2}]+...+B[1,j_{k}]=1.0;$$

• •

  thus if in step 1.1 above $Y_{\alpha}$ is the vector of approximations $y$ $./$ $\alpha$, then we have

  $$\mu_{1}[1]\leftarrow B[1,j_{1}]*Y_{\alpha}[j_{1}]+B[1,j_{2}]*Y_{\alpha}[j_{2}]+...+B[1,j_{k}]*Y_{\alpha}[j_{k}],$$

  which means that $\mu_{1}[1]$ is obtained as a weighted sum of approximations;

• •

  the same explanation goes for the steps 2.1 and 2.2 from each iteration: in 2.1, $\tilde{y}$ is the situation of the current iteration, and thus the vector $y$ $./$ $\tilde{y}$ is the vector of correction coefficients;

• •

  to see how $\mu_{1}[1]$ is computed in 2.2, observe that it is computed by multiplying the first line of the matrix diag($\mu_{1}$) $*$ $B$ with the vector of corrections $y$ $./$ $\tilde{y}$; thus if $Y_{\tilde{y}}$ is the vector of corrections $y$ $./$ $\tilde{y}$, then

  $$\mu_{1}[1]\leftarrow B[1,j_{1}]*\mu_{1}[1]*Y_{\tilde{y}}[j_{1}]+B[1,j_{2}]*\mu_{1}[1]*Y_{\tilde{y}}[j_{2}]+...+B[1,j_{k}]*\mu_{1}[1]*Y_{\tilde{y}}[j_{k}],$$

  which means that at each iteration $\mu_{1}[1]$ is obtained as a weighted sum of approximations (the weights being again $B[1,j_{1}]$, $B[1,j_{2}]$,…, $B[1,j_{k}]$).

Let $n_{x}=n_{y}=n_{d}=n_{v}=2$ (so $N=M=4$) and

$$\mu_{1}=\left[\begin{array}[]{cc}1.0&3.0\\ 2.0&4.0\end{array}\right]$$

the matrix to be reconstructed, written in vector format as

$$\mu_{1}=\left[\begin{array}[]{c}1.0\\ 2.0\\ 3.0\\ 4.0\end{array}\right].$$

We consider that the sinogram $y$ is obtained without any noise using the views V1 and V2, the view V1 at 0 degrees and the view V2 at a rotation around the object of 90 degrees, as in Figure 1; if each of the 4 squares of $\mu_{1}$ is of size 1.0 by 1.0, and the distance from X-ray source S to center of rotation (i.e. the center of reconstruction $\mu_{1}$) is 2.0, and the distance from center of rotation to the detectors’ line is also 2.0, and the detectors’ line is of length 4.0, this means that:

• •

  in view V1, $\mu_{1}[1]$ and $\mu_{1}[2]$ are each traversed by a full area of 1.0, whereas $\mu_{1}[3]$ and $\mu_{1}[4]$ are each traversed by an area of 0.75;

• •

  in view V2, $\mu_{1}[2]$ and $\mu_{1}[4]$ are each traversed by a full area of 1.0, whereas $\mu_{1}[1]$ and $\mu_{1}[3]$ are each traversed by an area of 0.75.

Thus

$$A=\left[\begin{array}[]{cccc}1.0&0.0&0.75&0.0\\ 0.0&1.0&0.0&0.75\\ 0.75&1.0&0.0&0.0\\ 0.0&0.0&0.75&1.0\end{array}\right],y=\left[\begin{array}[]{c}3.25\\ 5.0\\ 2.75\\ 6.25\end{array}\right],$$

$$\alpha=\left[\begin{array}[]{c}1.75\\ 1.75\\ 1.75\\ 1.75\end{array}\right],\beta=\left[\begin{array}[]{c}1.75\\ 2.0\\ 1.5\\ 1.75\end{array}\right],$$

$$B=\left[\begin{array}[]{cccc}\frac{1.0}{1.75}&0.0&\frac{0.75}{1.75}&0.0\\ 0.0&\frac{1.0}{2.0}&\frac{1.0}{2.0}&0.0\\ \frac{0.75}{1.5}&0.0&0.0&\frac{0.75}{1.5}\\ 0.0&\frac{0.75}{1.75}&0.0&\frac{1.0}{1.75}\end{array}\right];$$

the first line of $A$ corresponds to D1,V1, the second line to D2,V1, the third line to D1,V2, and the fourth line to D2,V2; also note that the fact that all values in $\alpha$ are equal is particular of this example (in other geometries they could be all different).

The SbIR algorithm works as follows (look only at $\mu_{1}[1]$):

• •

  in 1.1, $\mu_{1}[1]$ is initialized by multiplying the first line of $B$ with the vector $y$ $./$ $\alpha$, so

  $$\mu_{1}[1]=\frac{1.0}{1.75}*\frac{3.25}{1.75}+\frac{0.75}{1.75}*\frac{2.75}{1.75}=1.734;$$

  note that $\frac{3.25}{1.75}$ is the average sinogram measurement per unit of area in the D1,V1 pair and $\frac{2.75}{1.75}$ is the average sinogram measurement per unit of area in the D1,V2 pair; also note that $\mu_{1}[1]$ is traversed in all Detector,View pairs for a total area of 1.75 (this is $\beta[1]$), with an area of 1.0 traversed in the D1,V1 pair and an area of 0.75 in the D1,V2 pair;

• •

  in 2.1, in the first iteration, the situation is

  $$\tilde{y}=\left[\begin{array}[]{c}3.769\\ 4.662\\ 3.514\\ 5.3\end{array}\right].$$

• •

  in 2.2, in the first iteration, $\mu_{1}[1]$ is obtained by multiplying the first line of diag($\mu_{1}$) $*$ $B$ with the vector $Y_{\tilde{y}}$:

  $$\mu_{1}[1]=\frac{1.0}{1.75}*1.734*\frac{3.25}{3.769}+\frac{0.75}{1.75}*1.734*\frac{2.75}{3.514}=1.434;$$

  before this first iteration the situation in $\mu_{1}$ corresponding to the D1,V1 pair is 3.769, this means the elements that are traversed in D1,V1 need to be corrected by multiplication with $\frac{3.25}{3.769}$, so this is where the first approximation $1.734*\frac{3.25}{3.769}$ of $\mu_{1}[1]$ comes from; before this first iteration the situation in $\mu_{1}$ corresponding to the D1,V2 pair is 3.514, this means the elements that are traversed in D1,V2 need to be corrected by multiplication with $\frac{2.75}{3.514}$, so this is where the second approximation $1.734*\frac{2.75}{3.514}$ of $\mu_{1}[1]$ comes from;

• •

  as we can observe, $\mu_{1}[1]$ converges to its real value 1.0;

In the first example we consider the NCAT phantom studied also in other works [6] and shown in Fig. 2 (a) (this phantom simulates a thoracic cross-section); in this example we add noise to the sinogram by considering $rnl=0.005$ for both SbIR and the GPBB, UPN, UPN0 algorithms. Also the number of iterations is 512 for all algorithms (SbIR, GPBB, UPN, UPN0).

The reconstruction with the SbIR algorithm is shown in Fig. 2 (c), where in Fig. 2 (b) is shown the reconstruction obtained after the initialization step is applied. The reconstructions with the GPBB, UPN, UPN0 algorithms are shown in Fig. 3 when the regularization parameter $\alpha$ is 0.5, 5.0, or 50.0. In this example, the best reconstructions obtained with the GPBB, UPN, UPN0 algorithms are obtained when $\alpha=5.0$ (this means the middle row (b), (e) and (h) in Fig. 3).

We calculated the SSIM (Structural Similarity Index), PSNR (Peak Signal-to-Noise Ratio), MSE (Mean Square Error) for all the reconstructions: high values of SSIM and PSNR means higher accuracy while high values of MSE means lower accuracy; for the SbIR reconstruction they are SSIM:0.9018, PSNR: 25.6149, MSE: 0.0027; for GPBB, UPN, UPN0 they are listed in Table 1; as we can see the only case when the reconstruction with the SbIR algorithm is better is the one shown in Fig. 3 (c).

A problem with the GPBB, UPN, UPN0 algorithms is that for a low value of $\alpha$ the reconstruction is of low quality, whereas for a high value of $\alpha$ the reconstruction tends to become ”foggy” and with visible patches; in example 1, the value $\alpha=5.0$ seems to be a good choice.

We applied the same experimental setup but this time with $rnl=0.001$ (so with lower level of noise) to the abdominal phantom shown in Fig. 4, proposed by the Institute of Medical Physics, Friedrich-Alexander-University Erlangen-Nurnberg, Erlangen [23]; the reconstruction with UPN for this phantom was also best for $\alpha=5.0$ as in the case of example 1, whereas GPBB and UPN0 were best for $\alpha=0.5$ in this case. The conclusion would be that $\alpha=5.0$ works best in general for the UPN algorithm, for reconstructions of size 512 by 512 pixels; for GPBB and UPN0, the choice of $\alpha$ is more unclear.

In this second example we show that setting $\alpha=5.0$ could lead to reconstruction obtained by the UPN algorithm of lower quality than the one obtained with the SbIR algorithm. Consider the image shown in Fig. 5 (a), which corresponds to a cross-section from a pipe welding process, image acquired during 2015 at the Institute of Non-destructive Testing, Tomsk Polytechnic University, Tomsk (this is an image from a non-destructive material testing experiment instead of medical imaging, but it is equally important since reducing the number of views yields an overall reduction of the experiment time); in this example, we add noise to the sinogram by considering $rnl=0.001$ for both SbIR and the GPBB, UPN, UPN0 algorithms. Also the number of iterations is 512 for all algorithms (SbIR, GPBB, UPN, UPN0).

The reconstruction with the SbIR algorithm is shown in Fig. 5 (c), where in Fig. 5 (b) is shown the reconstruction obtained after the initialization step is applied. The reconstructions with the GPBB, UPN, UPN0 algorithms are shown in Fig. 6 when the regularization parameter $\alpha$ is 0.5, 5.0, or 50.0. The best reconstructions obtained with the GPBB, UPN, UPN0 algorithms are obtained when $\alpha=0.5$ (this means the first row (a), (d) and (g) in Fig. 6).

We calculated the SSIM, PSNR, MSE for all the reconstructions; for the SbIR reconstruction they are SSIM: 0.8820, PSNR: 35.5374, MSE: 2.7942e-04; for GPBB, UPN, UPN0 they are listed in Table 2; as we can see the reconstruction with the SbIR algorithm is better than all 6 reconstructions from Fig. 6 (b), (c), (e), (f), (h), (i), but worse than the reconstructions corresponding to $\alpha=0.5$ shown in Fig. 6 (a), (d), (g). Thus, this second example shows that even if in general the TV regularization techniques converge faster than the SbIR algorithm, there are cases as the one shown here where the SbIR algorithm provides a better reconstruction. Also, the SbIR algorithm has no parameters to be set as in the case of the TV regularization techniques.

The reconstructions have been obtained on a conventional PC with an i7-6500u @ 2.5 GHz 4-cores processor, and 8 GB of memory; only one core out of the four available has been used for each of the algorithms SbIR, GPBB, UPN, UPN0. The SbIR algorithm was implemented in C language in Visual C++ 2015; the GPBB, UPN, UPN0 algorithms of the TVReg package are implemented in C language. For the reconstructions shown in the two examples (all with 512 iterations) the computation time was around 150 seconds for the SbIR algorithm and between around 100 and 250 seconds for the GPBB, UPN, UPN0 algorithms. In [22], the complete code of our implementation in Visual C++ 2015 (using Win32 application, without any MFC or ATL code) of the SbIR algorithm can be found. In the code:

1. 1.

   ID_X_ITERATIVEMETHOD1 is the id of the ’Iterative Method 1’ menu option that runs the SbIR algorithm, with the geometry defined by the global variables.

The convergence of the SbIR algorithm, after a number of studies we have done, seems to behave as follows: first the SbIR algorithm spends a number (that depends on the specific input and noise in the sinogram) of iterations to stabilize within a certain interval, after which the rest of the iterations are spent to converge to the optimal solution. For the second example (the pipe welding process image) that is run with 512 iterations, the convergence behavior is shown in Fig. 7: the red graph shows the sum

$$\sum_{i=1}^{M}y[i]$$

of all sinogram values, which for the second example it is $8.987062\times 10^{6}$; and in blue graph it is plotted the sum

$$\sum_{i=1}^{M}\tilde{y}[i]$$

of all situation values, at each of the 512 iterations. We observe that in this plot the SbIR algorithm spends around 100 iterations to stabilize within a certain interval, and after that starts to converge to an approximate solution.

We have tried other initial solutions for the SbIR algorithm (like the solution given by the Filtered Back-Projection algorithm) to see if with other initializations it converges faster, but it seems that its initialization as shown in this work works best.