Equation (5.62) can be written in a matrix form as

$y=Ax\left(5.63\right)$

where y is the measurement vector, x is the image vector, and A is a M × N projection matrix.

5.5Combined Reconstruction

Combining the methods presented in the above three sections, some new techniques can be formulated. The combination can happen in the formula derivation, method implementation, or in applications. Some examples of combined reconstruction include:

(1)Iterative transform.

(2)Iterative reconstruction re-projection.

(3)Reconstruction using angular harmonics.

(4)Reconstruction by polynomial fit to projection.

In this section, the iterative transform method (also called continuous ART) is presented. The way in which the algorithm for discrete data is derived from the reconstruction formula qualifies it as a transform method, but its iterative nature and its representation of the image make the algorithm similar to the iterative algorithms derived by ART (Lewitt, 1983).

From Figure 5.5, suppose f(x, y) is 0 outside of Q and L(s, θn) is the length of the intersect segment of line (s, θn) with Q. The task of reconstruction can be described as reconstructing a function f(x, y), given its projections g(s, θn) for all real numbers s and a set of N discrete angles θn. The reconstruction formula operates on a “semi-discrete” function of g(s, θn) containing two variables, where θn is a discrete variable and s is a continuous variable.

For i ≥ 0, the (i + 1)-th step of the algorithm produces an image f(i + 1)(x, y) from the present estimation f(i)(x, y) given by

$f^{\left(i+1\right)}\left(x,y\right)=\{\begin{array}{c}0\text{if}f\left(x,y\right)\cancel{\in }\text{Q}\\ f^{\left(i\right)}\left(x,y\right)+\frac{g\left(s,{\theta }_n\right)-g^{\left(i\right)}\left(s,{\theta }_n\right)}{L\left(s,{\theta }_n\right)}\end{array}\left(5.66\right)$

where n = (i mod N) + 1, s′ = x cos θn + y sin θn, g(i)(s, θn) is the projection of f(i)(x, y) along the direction θn, and f(0)(x, y) is the given initial function. The sequence of images f(i)(x, y) generated by such a procedure is known to converge to an image which satisfies all the projections.

The physical interpretation of eq. (5.66) is as follows. For a particular angle θn, find the projection at angle θn of the current estimation of the image, and subtract this function from the given projection at angle θn to form the residual projection as a function of s for this angle. For each s of the residual, apply a uniform adjustment to all image points that lie on the line (s, θn), where the adjustment is such that g(i + 1)(s, θn) agrees with the given g(s, θn). Applying such an adjustment for all s results in a new image, which satisfies the given projection for this angle. The process is then repeated for the next angle in the sequence. It is often useful to view each step as a one-view back-projection of a scaled version of the residual.

Equation (5.66) specifies operations on “semi-discrete” functions. For a practical computation, an algorithm that operates on discrete data is required. Following the method used in Section 5.2, to estimate g(s, θn) from its samples g(mΔs, θn), an interpolating function q(·) containing one variable is introduced

$g\left(s,{\theta }_n\right)\approx {\displaystyle \sum _{m-M^{-}}^{M^{+}}g\left(m\text{Δs,}{\theta }_n\right)q\left(s-m\text{Δs}\right)}\left(5.67\right)$

Similarly, to estimate f(x, y) from its samples f(kΔx, lΔy), a basis function B(x, y) can be introduced that acts as an interpolating function of the two variables

$f\left(x,y\right)\approx {\displaystyle \sum _{k-K}^{K^{+}}{\displaystyle \sum _{l-L^{-}}^{L^{+}}f\left(k\text{Δ}x,l\text{Δ}y\right)B\left(x-\text{Δ}x,y-l\text{Δ}y\right)}}\left(5.68\right)$

Replacing f(x, y) by f(i)(x, y) in eq. (5.68) and substituting this expression in eq. (5.15), it gets

$g^{\left(i\right)}\left(s,\theta \right)={\displaystyle \sum _{k,l}^{}{f}^{\left(i\right)}\left(k\Delta x,l\Delta y\right)G_{k,l}^{\left(B\right)}\left(s,\theta \right)}\left(5.69\right)$

where

$G_{k,l}^{\left(B\right)}\left(s,\theta \right)={\displaystyle \int B\left(s\times \cos \theta -t\times \sin \theta -k\text{Δ}x,s\times \sin \theta +t\times \cos \theta -l\text{Δ}y\right)\text{d}t}\left(5.70\right)$

According to the principles of the interpolation stated above, the continuous ART can be discretized. Deriving from eq. (5.67) and eq. (5.69), eq. (5.66) can be approximated as

$\begin{array}{cc}{f}_{k,l}^{\left(i+1\right)}=\{\begin{array}{cc}0& if\left(k\text{Δ}x,l\text{Δ}y\right)\cancel{\in }\text{Q}\\ f_{k,l}^{\left(i\right)}& \frac{{\displaystyle \sum _m^{}\left[g\left(m\text{Δ}s,{\theta }_n\right)-{\displaystyle \sum _{k,l}^{}{f}_{k,l}^{\left(i\right)}\times G_{k,l}^{\left(B\right)}\left(m\text{Δ}s,{\theta }_n\right)}\right]q\left[s_{k,l}\left({\theta }_n\right)\right]-m\text{Δ}s}}{L\left[s_{k,l}\left({\theta }_n\right),{\theta }_n\right]}\end{array}& \left(5.71\right)\end{array}$

where n = (i mod N) + 1, and

$f_{k,l}^{\left(i\right)}=f^{\left(i\right)}\left(k\text{Δ}x,l\text{Δ}y\right)\left(5.72\right)$

5.6Problems and Questions

Figure Problem 5-1

5-2A 5 × 5 matrix is composed of elements computed from 99/[1+(5|x|)1/2+(9|y|)1/2], which are rounded to the nearest integer. The values of x and y are from –2 to 2. Compute the projections along the four directions with θ = 0°, 30°, 60°, and 90°.

5-3Reconstructing a density distribution f(x, y) from three projections is certainly not enough. Some people believe that if a priori knowledge f(x, y) = ∑δ(x – x′, y – y′) were available, the exact reconstruction from just three projections would be possible. However, this rather restrictive condition implies that all point masses are equal and none of them are negative. Others say, “If you can go that far, can’t you handle positive impulse patterns of unequal strength?” What is your opinion?

5-4The gray levels of a 4 × 4 image are shown in Figure Problem 5-4, four of which are known (the question mark means unknown). In addition, its projection on the X axis is P(x) = [? ? 15 5]T, and on the Y axis it is P(y) = [9 ? 27 ?]T. Can you reconstruct the original gray-level image?

$\begin{array}{cccc}?& 3& ?& ?\\ 0& ?& ?& ?\\ 6& ?& ?& ?\\ ?& 3& ?& ?\end{array}$

Figure Problem 5-4

Figure Problem 5.5

5-7*Prove the validity of the projection theorem for Fourier transforms.

5-8Use the projection data obtained in Problem 5-1 to reconstruct the original images with reconstruction from Fourier inversion.

5-9What are the relations and differences among the methods of reconstruction based on back projection?

5-11Implement the ART with relaxation technique in Section 5.4 and verify the data in Figure 5.20.

5-12Perform a reconstruction according to the projection data from Problem 5-1 using ART without relaxation and ART with relaxation, respectively. Compare the results with the original images in Problem 5-1.

5.7Further Reading

1.Modes and Principles

–In history, the first clinical machine for detection of head tumors, based on CT, was installed in 1971 at the Atkinson Morley’s Hospital, Wimbledon, Great Britain (Bertero and Boccacci, 1998). The announcement of this machine, by G.H. Hounsfield at the 1972 British Institute of Radiology annual conference, is considered the greatest achievement in radiology since the discovery of X-rays in 1895.

–In 1979, G.H. Hounsfield shared with A. Cormack the Nobel Prize for Physiology and Medicine (Bertero and Boccacci, 1998). The Nobel Prizes for Chemistry of 1981 and 1991 are also related to it (Herman, 1983), (Bertero and Boccacci, 1998), and (Committee, 1996).

–The phenomenon of nuclear resonance was discovered in 1946 by F. Bloch and M. Purcell, respectively. They shared the Nobel Prize for Physics in 1952 (Committee, 1996). The first MRI image was obtained in 1973.

–Complementary reading on the history and development of image reconstruction from projections can be found in Kak and Slaney (2001).

2.Reconstruction by Fourier Inversion

–An introduction and some examples of reconstruction by Fourier inversion can be found in Russ (2002).

–Another reconstruction method using 2-D Fourier transform is rho-filtered layergrams Herman (1980).

–More discussion and applications of phantoms can be found in Moretti et al. (2000).

–Another head phantom for testing CT reconstruction algorithms can be found in Herman (1980).

3.Convolution and Back-Projection

–The advantages of convolution and back-projection include that they are more suitable than other methods for data acquisition using a parallel mode (Herman, 1980).

–A recent discussion on back-projection algorithms and filtered back-projection algorithms can be found in Wei et al. (2005).

–More mathematical expressions on cone-beam geometry and 3-D image reconstruction can be found in Zeng (2009).

4.Algebraic Reconstruction

–The methods for algebraic reconstruction have many variations and further details can be found in Herman (1980), Censor (1983), Lewitt (1983), Zeng (2009), and Russ (2016).

5.Combined Reconstruction

–More combined reconstruction techniques, particularly second-order optimization for reconstruction and the non-iterative finite series-expansion reconstruction methods can be found in Herman (1980, 1983).