SIAM News Blog
SIAM News
Print

What Can You See with Computed Tomography?

By Per Christian Hansen, Jakob Sauer Jørgensen, William R.B. Lionheart, and Eric Todd Quinto

Computed Tomography: Algorithms, Insight, and Just Enough Theory—edited by Per Christian Hansen, Jakob Sauer Jørgensen, and William R.B. Lionheart, with contributions from Martin S. Andersen, K. Joost Batenburg, Yiqiu Dong, Eric Todd Quinto, and Jan Sijbers—was published by SIAM in 2021. The book, which is part of the Fundamentals of Algorithms series, illustrates how theoretical insight and numerical computation go hand in hand in the analysis of algorithms for computed tomography (CT).

The following text focuses on CT problems with limited data and highlights a central theme in inverse problems: What can we reliably expect to recover from specific measurements? The text is based on excerpts from chapters 5, 8, and 10 (where more illustrations can be found) and has been modified slightly for clarity.

We begin by introducing a few basic elements of CT. Without loss of generality, we assume that we are scanning an object that is contained in a unit disk of radius \(1\). The object’s attenuation coefficient is characterized by the nonnegative function \(f(x_1, x_2)\). Let \(L_{\theta,s}\) denote a line at angle \(\theta\) and a distance \(s\) from the origin. We define the projection that expresses all line integrals along \(L_{\theta,s}\) at fixed angle \(\theta\) as

\[g_\theta(s)=\int_{L_{\theta,s}}f(x_1,x_2)d\ell \qquad \textrm{for} \qquad s\in[-1,1].\tag1\]

The restriction of \(s\) to \([-1,1]\) is possible because all lines for \(s\) that fall outside of this interval do not intersect with the object, meaning that their line integrals are simply equal to \(0\). We then define the Radon transform of \(f\) as the projections at all angles, i.e.,

\[\mathcal{R}[f](\theta,s)=g_\theta(s) \qquad \rm{for} \qquad \theta \in [0^{\circ}, 360^{\circ}[\tag2\]

and \(s\in[-1,1]\). The output of the Radon transform is also known as the sinogram due to distinct sinusoidal bands that appear when one views it as an image.

Complete data are taken over lines for all angles \(\theta\in[0^{\circ},360^{\circ}[\) and displacements \(s\in[-1,1]\), while limited data are taken over smaller sets of lines. Here we mention three common types of limited data:

  • Limited-angle data arise when the angles \(\theta\) are restricted. For example, limited-angle data would ensue if \(\theta\) was restricted to \([90^\circ, 180^\circ]\cup[270^\circ, 360^\circ]\) but all values of \(s\in[-1,1]\) were allowed.
  • Region-of-interest data arise when lines only pass through a specific region of interest. For instance, if the region of interest is a ball of radius \(1/2\) at the origin, then the corresponding region-of-interest data would be over lines for \(\theta\in[0^\circ,360^\circ[\) and \(s\in[-1/2,1/2]\).
  • Exterior data arise when we limit the data to avoid a specific region. If the excluded region is the ball of radius \(1/2\) that is centered at the origin, for example, then the exterior data would exist over all lines that do not meet this ball, i.e., lines for \(\theta\in[0^\circ,360^\circ[\) and \(s\in[-1,-1/2]\cup[1/2,1]\).

Figure 1. Illustration of limited-angle reconstruction in Example 1. 1a. A disk of radius \(1/2\) that is centered at the origin. 1b. A limited-angle sinogram. 1c. The corresponding reconstruction.

Example 1 

Now we illustrate a simple limited-angle reconstruction. Figure 1 depicts a disk of radius \(1/2\) that is centered at the origin, a limited-angle sinogram for \(\theta\in[90^\circ,180^\circ]\cup[270^\circ,360^\circ]\), and a limited-angle reconstruction of the disk, computed as described in chapter 6 of Computed Tomography. Here, the limited-angle data domain solely consists of lines for \(\theta\in[90^\circ,180^\circ]\cup[270^\circ,360^\circ]\) and \(s\in[-1,1]\). The sinogram in Figure 1b is black in the missing range and only shows the data values in the limited data domain. The magenta vertical lines at \(\theta=90^\circ\), \(\theta=180^\circ\), \(\theta=270^\circ\), and \(\theta=360^\circ\) represent the data domain’s outer boundaries.

The reconstruction in Figure 1c should look like the disk in Figure 1a, but it appears markedly different. The reconstruction does correctly contain circle boundaries with positive slopes in the second and fourth quadrants. The lines for \(\theta\) between \(90^\circ\) and \(180^\circ\) have positive slopes, so the well-reconstructed object boundaries are tangent to lines in the data domain; these are the visible boundaries. However, the negatively sloped boundaries in the first and third quadrants are not evident in the reconstruction. Because of their negative slopes, these boundaries are tangent to no lines in the data domain and are thus invisible.

For \(s=\pm 1/2\), we see vertical artifacts along the lines with \(\theta=0^\circ\) and horizontal artifacts along the lines with \(\theta=90^\circ\). Recall that the data’s outer boundary consists of lines with \(\theta=90^\circ\), \(180^\circ\), \(270^\circ\), and \(360^\circ\) (see Figure 1b). Therefore, the artifact lines are tangent to the object and exist in the outer boundary of the data domain. \(\blacksquare\)

Figure 2. The first 20 right singular vectors \(\boldsymbol{v}_i\) of the matrix \(\boldsymbol{A}\) for the limited-angle problem.
Using the discretization methods from chapter 9 of Computed Tomography, we obtain a linear algebra formulation of the CT reconstruction problem that takes the form \(\boldsymbol{Ax}=\boldsymbol{b}\). In this form, vectors \(\boldsymbol{x}\) and \(\boldsymbol{b}\) respectively represent the image and sinogram. Let matrix \(\boldsymbol{A}\) have dimensions \(m \times n\) and assume that \(m \geq n\) for the sake of simplicity. The singular value decomposition (SVD) then takes the form

\[\boldsymbol{A}=\sum^n_{i=1}\boldsymbol{u}_i\sigma_i\boldsymbol{v}_i^T,\tag3\]

where \(\sigma_i\) are the singular values and \(\boldsymbol{u}_i\) and \(\boldsymbol{v}_i\) are the corresponding left and right singular vectors. Further details about the SVD are available in chapter 10.

Example 2 

Now we have computed the leading singular values and vectors of matrix \(\boldsymbol{A}\) from Example 1. Figure 2 illustrates some of the right singular vectors \(\boldsymbol{v}_i\), which have distinctive features that are characteristic of this specific limited-angle problem. In particular, all of the reshaped vectors are dominated by features that are elongated at a \(45^\circ\) angle with respect to the horizontal (the features for the complete-data problem have no dominating direction).

At this point, we know from the analysis in chapter 8 of Computed Tomography that we can only expect to recover structures and edges along the projection angles; other features are much harder to reconstruct without a priori information. Our singular vectors are completely in accordance with this fact. Since the reconstruction is a weighted sum over the right singular vectors \(\boldsymbol{v}_i\), it will inherit the dominating features of these vectors — in this case, features and structures at roughly \(45^\circ\) angles. \(\blacksquare\)

Due to ill conditioning of \(\boldsymbol{A}\), high-frequency components that relate to the smallest singular values tend to contaminate naïve CT reconstructions. A simple approach that reduces this noise’s influence involves merely discarding the SVD coefficients that correspond to the small singular values. Doing so leads to the truncated SVD (TSVD) solution

\[x_k=\sum^k_{i=1}\frac{\boldsymbol{u}^T_i\boldsymbol{b}}{\sigma_i}\boldsymbol{v}_i, \qquad k<n.\tag4\]

Example 3 

We finish by computing TSVD solutions to the limited-angle problem from Example 2 (see Figure 3). As \(k\) increases, contributions from higher-frequency singular vectors \(\boldsymbol{v}_i\) are included, producing an increasingly sharper image. As Figure 2 shows, certain features dominate the right singular vectors \(\boldsymbol{v}_i\); as a result, the reconstructions will also account for these features. We can thus indeed reconstruct the features and edges of the exact image at angles that are around \(45^\circ\). On the other hand, features and edges at approximately \(-45^\circ\) angles are almost completely absent from the reconstructions. \(\blacksquare\) 

Figure 3. Truncated singular value decomposition (TSVD) solutions to the limited-angle problem for select values of the truncation parameter \(k\), as well as the exact solution.

While there is no straightforward way to improve TSVD, more advanced reconstruction methods come to the rescue. These methods include algebraic iterative reconstruction methods (in chapter 11) and regularization techniques to incorporate prior information that compensates for incomplete data (in chapter 12). The latter rely on efficient optimization methods that are described in chapter 13 of Computed Tomography.


Enjoy this passage? Visit the SIAM Bookstore to learn more about Computed Tomography: Algorithms, Insight, and Just Enough Theory and browse other SIAM titles.

Per Christian Hansen is a professor of scientific computing at the Technical University of Denmark (DTU) and a SIAM Fellow with a research track in numerical methods for inverse problems. Jakob Sauer Jørgensen is a senior researcher at DTU and a Presidential Fellow at the University of Manchester. His work focuses on algorithms and software for computed tomography and uncertainty quantification. William R.B. Lionheart is a professor of applied mathematics at the University of Manchester who specializes in applied and theoretical inverse imaging problems in medicine, materials science, security, and radar. Eric Todd Quinto is the Robinson Professor of Mathematics at Tufts University. He uses microlocal analysis to understand strengths and weaknesses that are inherent to limited data problems in a range of tomographic modalities.

blog comments powered by Disqus