Machine Learning for Tomographic Imaging. Professor Ge Wang

characteristics. The relaxation time has two main aspects: T1 (longitudinal relaxation time) and T2 (transverse relaxation time). While CT images enjoy excellent bone–tissue and air–tissue contrasts, MRI images give rich soft tissue contrasts and functional information.

      In addition to x-ray CT and MRI, there are multiple other tomographic imaging modalities. Importantly, nuclear emission imaging concentrates on physiological functions such as blood flow and metabolism. In this category, positron-emission tomography (PET) and single-photon emission computed tomography (SPECT) are the primary modes of nuclear tomography (Leigh et al 2002, Zeng 2009).

      A PET system detects pairs of gamma rays generated by a positron-emitting radionuclide, such as fluorine-18, which is introduced into a patient inside a biologically active molecule called a radioactive tracer. A pair of gamma ray photons is emitted in opposite directions from the radioactive tracer. Hence, the PET scanner utilizes their linearity and simultaneity to detect the emission event along a corresponding line of response defined by two detectors facing each other. Specifically, if two gamma ray photons are simultaneously (within a time window of 10–20 ns) detected, then a signal is recorded for the line of response. From these data, a distribution of the radioactive tracer can be reconstructed. Similarly, SPECT reconstructs distributions of radioactive tracers from gamma ray induced signals, but these tracers only emit gamma ray photons individually, and a metallic collimator is needed to determine the line path along which a gamma ray photon travels. In both the PET and SPECT cases, the tracer nuclides require a patient to contain a radiation source in his/her body, and the detector measures radiation-induced data externally for image reconstruction.

      Furthermore, we can also utilize ultrasound and light waves for imaging purposes, which are called ultrasound imaging and optical imaging, respectively. Generally speaking, for any physical measurement on an object of interest, as long as the data do not directly reflect structural features, a so-called ‘inverse’ process will be needed to estimate these features from the indirect data. Tomography is a very important class of inverse problems, targeting cross-sectional/volumetric image formation. In the next section, we will heuristically explain the CT problem as a specific example.

1.1.2 Radon transform and non-ideality in data acquisition

      The Radon transform describes the relationship between an underlying function and a simple indirect measurement process. First studied in 1917, Johann Radon sought to invert this transform and reconstruct the underlying function from these measurements. It is mathematically natural, as it applies an integral operator over a sub-space (for example, a line or a plane) in a high-dimensional space (for example, a 3D Euclidean space). Also, it is physically relevant, as an x-ray signal attenuated by an object can be easily converted into a line or planar integrals after practical approximations.

      Without loss of generality, let us consider a function defined on a plane and perform the Radon transform along lines in the plane. Then, the value of the Radon transform along an arbitrary line is equal to the following line integral:

      Rf(L)=∫Lf(x)dx,(1.1)

      where Rf denotes projection data that can be acquired by a CT scan and f represents an underlying function depicting linear attenuation coefficients inside the object. Technically, x-ray photons go through the object, are attenuated, and are then recorded by detector elements along a 1D array from various viewing angles. X-ray data can be put in a 2D array with respect to the viewing angle and the detector location, forming a so-called sinogram, as shown in figure 1.1.

      Figure 1.1. Radon transform as a simple example of indirectly measured data, from which an underlying function needs to be estimated or reconstructed.

      In this 2D case, the most commonly used analytical formula to recover f from its Radon transform is the filtered backprojection (FBP) formula. In the 3D case, the Radon transform typically gives planar integrals, instead of line integrals, and can be inverted using more complicated formulas. In a nutshell, the Radon transform is closely related to the Fourier transform. Thus, if all Radon data is available, then the inverse Radon transform is essentially the inverse Fourier transform. Alternatively, we can view each line or planar integral as a linear equation. With all available x-ray data, we have a system of linear equations. In principle, we can solve this linear system in some way to uncover all unknown pixel/voxel values, i.e. to reconstruct the underlying function/image. We will discuss specific technical details in chapter 4, which is dedicated to CT basics.

      In the ideal case, data indirectly measured by many tomographic imaging systems can often be approximated as linear combinations of unknown variables assuming neither noise nor bias. However, in reality there are many practical factors that prevent us from obtaining idealized data. During a measurement process, interactions between a physical probing method and an object to be reconstructed are often stochastic processes, and both the object and the imaging system can be time-variant, introducing uncertainties and inconsistencies in the data. For example, the measurement of x-ray or gamma ray photons contains inherent Poisson noise. Also, current-integrating x-ray detectors have electronic noise.

      Given the radiation risk of x-ray radiation, which might carry genetic, cancerous, and other risks, CT scanning with a reduced radiation dose has been a hot topic over the past decade. Currently, the medical community is striving for high-quality CT images at a lowest possible dose level, for example, in the Image Gently campaign on pediatric patients and the Image Wisely campaign on adult patients the ‘as low as reasonably achievable’ (ALARA) principle has been widely accepted. These efforts bring up the well-known low-dose CT challenge. The low-dose condition will severely degrade the image quality, since x-ray imaging is a quantum accumulation process.

      As another example, an MRI scan takes much longer than a CT scan does, and needs to be accelerated with fast MRI techniques. As a result, measured MRI data, known as Fourier space (or k-space) samples, cannot fully cover the Fourier space. Simply applying the inverse Fourier transform to the partially measured Fourier spectrum produces an image with strong artifacts, in particular when these data are further compromised with patient and/or physiological motion.

      In the above CT and MRI issues, which are among many imaging problems, analytic image reconstruction methods are not suitable, as they demand ideal data and are almost exclusively based on the Fourier formulation. To alleviate these types of problems, iterative image reconstruction methods are advantageous. Iterative reconstruction algorithms use optimization techniques, easily incorporate prior knowledge on the imaging process and the image information content, and solve the imaging problem iteratively at an increased computational cost. Compared to analytic reconstruction algorithms, iterative algorithms reduce image noise and artifacts effectively.

1.1.3 Bayesian reconstruction

      As mentioned before, tomography is nothing other than image reconstruction from measurement data indirectly related to hidden structures. In a simple case, such as the Radon transform, data b is related to an underlying image x through a linear mapping A. This mapping is called the forward process:

b=Ax.(1.2)

      In the case of a CT scan, b denotes the data acquired by the scan, x represents the underlying CT image, and A is the system matrix specific to a CT scanner and the imaging protocol used for the CT scan. Actually, equation (1.2) is a discrete form of the Radon transform. Now, our inverse problem is to calculate the image x from data b, given the imaging system matrix A.

      The analytical solution to the inverse problem in the Fourier domain has a close-form expression but cannot handle image noise and artifacts well. Another straightforward way to resolve this problem is to compute