2.2. FBP

FBP algorithms are widely used in CT when many projections acquired at >360 are used to reconstruct cross-sectional images. The number of projections typically ranges from a few hundred to approximately 1000. The Fourier central slice theorem is fundamental to FBP theory. In two-dimensional (2D) CT imaging, a projection of an object corresponds to sampling that object along the direction perpendicular to the X-ray beam in the Fourier space [13]. For many projections, information about the object is well sampled, and the object can be restored by combining the information from all projections.

In three-dimensional (3D) cone-beam imaging, information about the object in Fourier space is related to the Radon transform of the object. The relationship between the Radon transform and cone-beam projections has been well studied, and solutions to the cone-beam reconstruction have been provided [14, 15]. The Feldkamp algorithm generally provides a high degree of precision for 3D reconstruction images when an exact type of algorithm is employed [16]. Therefore, this method has been adopted for the image reconstruction of 3D tomography and multidetector cone-beam CT. A number of improved 3D reconstruction methods have been derived from the Feldkamp algorithm.

We denoted the intensity of the incident X-rays as I<sup>0</sup> ¼ ð Þ δ; c; d and that of the X-rays that passed through the structure at the location ð Þ <sup>δ</sup>; <sup>c</sup>; <sup>d</sup> as <sup>I</sup> <sup>¼</sup> ð Þ <sup>δ</sup>; <sup>c</sup>; <sup>d</sup> . The image data dF <sup>¼</sup> ð Þ <sup>δ</sup>; <sup>c</sup>; <sup>d</sup> were calculated as follows:

$$d^{\mathbb{F}}(\delta, c, d) = \ln \frac{I\_0(\delta, c, d)}{I(\delta, c, d)} \tag{1}$$

The original projection data are used to generate an FBP image:

$$f^{\text{FDK}} = \int\_{-\delta}^{\delta} \frac{R^2}{\Omega^2} \hat{d}F(\delta, c, d)D\delta \tag{2}$$

$$\mathcal{U}(X, Y, \delta) = R + X\cos\delta + Y\sin\delta$$

$$c(X, Y, \delta) = R\frac{-X\sin\delta + Y\cos\delta}{R + X\cos\delta + Y\sin\delta}$$

$$d(X, Y, Z, \delta) = Z\frac{R}{R + X\cos\delta + Y\sin\delta}$$

$$\hat{d}F(\delta, c, d) = \left(\frac{R}{\sqrt{R^2 + c^2 + d^2}} d^\delta(\delta, c, d)\right) \* \mathcal{g}^P(\gamma)$$

where the projection angle <sup>δ</sup>, fan-angle <sup>γ</sup>, source trajectory <sup>R</sup>, acquisition angle <sup>θ</sup>, and gpð Þ <sup>γ</sup> represent a convolution with the Ramachandran–Lakshminarayanan filter. The Ramachandran– Lakshminarayanan (ramp) filter is shown below:

$$\mathbf{g}^p(\boldsymbol{\gamma}) = \left(\frac{\boldsymbol{\gamma}}{\sin \alpha}\right)^2 \tag{3}$$

#### 2.3. IR

Interest in digital tomosynthesis (DT) and its clinical applications has been revived by recent advances in digital X-ray detector technology. Conventional tomography technology provides planar information about an object from its projection images. In tomography, an X-ray tube and an X-ray film receptor lie on either side of the object. The relative motion of the tube and film is predetermined on the basis of the location of the in-focus plane [1]. A single image plane is generated by a scan, and multiple scans are required to provide a sufficient number of planes to cover the selected structure in the object. DT acquires only one set of discrete X-ray projections that can be used to reconstruct any plane of the object retrospectively [2]. The technique has been investigated in angiography and in the imaging of the chest, hand joints, lungs, teeth, and breasts [3–8]. Dobbins et al. [9] reviewed DT and showed that it outperformed planar imaging to a statistically significant extent. Various types of DT reconstruction methods

Current DT mainly involves image acquisition/reconstruction using polychromatic imaging. Material decomposition or virtual monochromatic image processing using dual energy (DE) has been studied in CT, and many basic and clinically useful applications have been reported [10–12]. Similar to CT, it can be expected that DT also benefits from image quality improvements. In this chapter, the fundamental image quality characteristics of various reconstruction algorithms (including a state-of-the-art reconstruction algorithm) using polychromatic imaging and virtual monochromatic DT imaging that were verified in phantom experiments are discussed.

Existing tomosynthesis algorithms can be divided into three categories: (1) backprojection

The BP algorithm is referred to as "shift-and-add" (SAA) [9], whereby projection images taken at different angles are electronically shifted and added to generate an image plane focused at a certain depth below the surface. The projection shift is adjusted such that the visibility of the features in the selected plane is enhanced, whereas those in other planes are blurred. By using a digital detector, the image planes at all depths can be retrospectively reconstructed from one set of projections. The SAA algorithm is valid only if the motion of the X-ray focal spot is

FBP algorithms are widely used in CT when many projections acquired at >360 are used to reconstruct cross-sectional images. The number of projections typically ranges from a few hundred to approximately 1000. The Fourier central slice theorem is fundamental to FBP theory. In two-dimensional (2D) CT imaging, a projection of an object corresponds to sampling that object along the direction perpendicular to the X-ray beam in the Fourier space [13]. For

(BP), (2) filtered backprojection (FBP), and (3) iterative reconstruction (IR) algorithms.

have been explored.

68 Medical Imaging and Image-Guided Interventions

2. DT reconstruction

parallel to the detector.

2.1. BP

2.2. FBP

An IR algorithm performs the reconstruction in a recursive fashion [17, 18], unlike the one-step operation in backprojection and FBP algorithms. During IR, a 3D object model is repeatedly updated until it converges to the solution that optimizes an objective function. The objective function defines the criteria of the reconstruction solution.

Algebraic reconstruction methods assume that the image is an array of unknowns, and the reconstruction problem can be established as a system of linear equations. The unknowns of this system are approximated with respect to the ray sums iteratively [17, 18]. In each iteration, current reconstruction is reprojected and updated according to how much it satisfies the projection measurements.

$$\mathbf{G}\_{\mathbf{x}} = \mathbf{f} \tag{4}$$

x !ð Þ<sup>k</sup> ¼ x !ð Þ <sup>k</sup>�<sup>1</sup>

mizes the probability of obtaining the measured projections.

<sup>C</sup>ð Þ<sup>k</sup> ðÞ¼ <sup>i</sup> <sup>1</sup> P<sup>n</sup>

<sup>x</sup>ð Þ <sup>k</sup>þ<sup>1</sup> ðÞ¼ <sup>i</sup> <sup>1</sup> P<sup>n</sup>

2.4. Total variation minimization reconstruction algorithm

<sup>j</sup>¼<sup>1</sup> <sup>α</sup>ð Þ <sup>i</sup>; <sup>j</sup>

the vector y(j)=[y(1), y(2), y(3),..., y(J)]:

i = 1, 2, …, n.

þ β

where Ω<sup>t</sup> is the set of indices of the rays sent from the tth projection direction.

1 Σ<sup>i</sup> <sup>∈</sup> <sup>Ω</sup><sup>t</sup> aij

The objective function in the maximum likelihood expectation–maximization (MLEM) algorithm is the likelihood function, which is the probability of obtaining the measured projections in a given object model. The solution of the MLEM algorithm is an object model that maxi-

The transition matrix αð Þ i; j represents the probability for an event generated in the area of the source covered by pixel i to be detected. The count registered by the detector is represented by

> X J

> > j¼1

<sup>j</sup>¼<sup>1</sup> <sup>α</sup>ð Þ <sup>i</sup>; <sup>j</sup> <sup>x</sup>ð Þ<sup>k</sup> ð Þ<sup>i</sup>

An iterative algorithm using total variation (TV)-based compressive sensing was recently developed for volume image reconstruction from a tomographic scan [21]. The image TV has been used as a penalty term in iterative image reconstruction algorithms [21]. The TV of an image is defined as the sum of the first-order derivative magnitudes for all pixels in the image. TV minimization is an image domain optimization method associated with compressed sensing theory [21]. As TV-minimization IR for image reconstruction, the TV-minimization IR technique is the SART [20] with algebraic IR for constraining the TV-minimization problem, which is called SART-TV [21]. In TV-minimization IR, adding a penalty to the data–fidelity–objective function tends to smooth out noise in the image while preserving edges within the image [21–25].

SART-TV minimizes the Rudin, Osher, and Fatemi (ROF) model [26], that is, SART-TV also takes into account the image information when minimizing the TV of the image. If only the TVminimization step was performed in the rest of the algorithms, the result would be a flat image; alternatively, the ROF model ensures that the image does not change considerably. The SART-TV optimal parameters include the iteration number for TV minimization [100 (np)] and the length of each gradient-descent step [50 (q)]. These optimal parameters have been

P<sup>n</sup>

X i ∈ Ω<sup>t</sup> aij bi � a ! i x !ð Þ <sup>k</sup>�<sup>1</sup>

P<sup>n</sup> <sup>j</sup>¼<sup>1</sup> aij

<sup>x</sup>ð Þ <sup>k</sup>þ<sup>1</sup> ðÞ¼ <sup>i</sup> <sup>x</sup>ð Þ<sup>k</sup> ð Þ<sup>i</sup> <sup>C</sup>ð Þ<sup>k</sup> ð Þ<sup>i</sup> (8)

y jð Þαð Þ i; j

<sup>i</sup>¼<sup>1</sup> <sup>α</sup>ð Þ <sup>i</sup>; <sup>j</sup> <sup>x</sup>ð Þ<sup>k</sup> ð Þ<sup>i</sup> <sup>α</sup>ð Þ <sup>i</sup>; <sup>j</sup> (9)

State-Of-The-Art X-Ray Digital Tomosynthesis Imaging http://dx.doi.org/10.5772/intechopen.81667

<sup>i</sup>¼<sup>1</sup> <sup>x</sup>ð Þ<sup>k</sup> <sup>α</sup>ð Þ <sup>i</sup>; <sup>j</sup> (10)

y jð Þ

P<sup>n</sup>

X J

j¼1

(7)

71

In the DT reconstruction, x represents the pixel values for x∈ ℜ<sup>N</sup>, whereas f represents the observed data for f ∈ ℜ<sup>M</sup>. The weighting value G ∈ ℜ<sup>M</sup>�<sup>N</sup> was created by using a forward projection, and the matrix G combines the submatrices of each projection. Gmn is defined as the length of intersection of the mth ray with the nth cell. DT reconstruction starts with an initial guess x<sup>0</sup> and computes x<sup>1</sup> by projecting x<sup>0</sup> onto the first hyperplane in Eq. 4. The algebraic reconstruction technique (ART) update procedure (or error correction) is shown below:

$$\overrightarrow{\mathbf{x}}^{(k)} = \overrightarrow{\mathbf{x}}^{(k-1)} + \beta \frac{b\_i - \overrightarrow{a}\_i \overrightarrow{\mathbf{x}}^{(k-1)}}{\left\| \overrightarrow{a}\_i \right\|^2} \overrightarrow{a}\_i \tag{5}$$

where ai !¼ ð Þ ai1; ai2;…; ain is the <sup>i</sup>th row of the projection matrix, and <sup>β</sup> is the relaxation parameter (1.0). a ! i � � � � � � 2 is the normalization factor. The update is performed for each projection measurement bi separately, that is, the kth iteration consists of a sweep through the m projection measurements. The algorithm iterates through the equations periodically; therefore, i = (i-1) mod (m) + 1.

The ART algorithm updates the image vector per ray such that the update satisfies only a single equation representing the corresponding ray integral. By contrast, the simultaneous IR technique (SIRT) algorithm updates the image vector after all equations are considered. The update procedure of SIRT is given in Eq. 6 according to [19].

$$\overrightarrow{\mathbf{x}}^{(k)} = \overrightarrow{\mathbf{x}}^{(k-1)} + \beta \frac{1}{\sum\_{i}^{m} a\_{i\bar{\jmath}}} \sum\_{i}^{m} a\_{i\bar{\jmath}} \frac{b\_{i} - \overrightarrow{a}\_{i} \overrightarrow{\mathbf{x}}^{(k-1)}}{\sum\_{j=1}^{n} a\_{i\bar{\jmath}}} \tag{6}$$

The simultaneous ART (SART) algorithm [20] is proposed as a combination of the ART and SIRT algorithms. SART updates the superior implementation of ART and is based on a simultaneous update of the current reconstruction similar to SIRT. In the SART algorithm, the update procedure is applied for all rays in a given scan direction (projection) instead of each ray separately (similar to conventional ART) or instead of all rays simultaneously (similar to SIRT). The SART update is expressed as follows:

State-Of-The-Art X-Ray Digital Tomosynthesis Imaging http://dx.doi.org/10.5772/intechopen.81667 71

$$\overrightarrow{\mathbf{x}}^{(k)} = \overrightarrow{\mathbf{x}}^{(k-1)} + \beta \frac{1}{\sum\_{i \in \Omega\_l} a\_{ij}} \sum\_{i \in \Omega\_l} a\_{ij} \frac{b\_i - \overrightarrow{a}\_i \overrightarrow{\mathbf{x}}^{(k-1)}}{\sum\_{j=1}^{n} a\_{ij}} \tag{7}$$

where Ω<sup>t</sup> is the set of indices of the rays sent from the tth projection direction.

The objective function in the maximum likelihood expectation–maximization (MLEM) algorithm is the likelihood function, which is the probability of obtaining the measured projections in a given object model. The solution of the MLEM algorithm is an object model that maximizes the probability of obtaining the measured projections.

The transition matrix αð Þ i; j represents the probability for an event generated in the area of the source covered by pixel i to be detected. The count registered by the detector is represented by the vector y(j)=[y(1), y(2), y(3),..., y(J)]:

$$\mathbf{x}^{(k+1)}(i) = \mathbf{x}^{(k)}(i)\mathbf{C}^{(k)}(i) \tag{8}$$

$$\mathbf{C}^{(k)}(i) = \frac{1}{\sum\_{j=1}^{n} \alpha(i,j)} \sum\_{j=1}^{l} \frac{y(j)}{\sum\_{i=1}^{n} \alpha(i,j)\mathbf{x}^{(k)}(i)} \alpha(i,j) \tag{9}$$

$$\mathbf{x}^{(k+1)}(\mathbf{i}) = \frac{1}{\sum\_{j=1}^{n} \alpha(\mathbf{i}, \mathbf{j})} \mathbf{x}^{(k)}(\mathbf{i}) \sum\_{j=1}^{l} \frac{y(\mathbf{j})\alpha(\mathbf{i}, \mathbf{j})}{\sum\_{i=1}^{n} \mathbf{x}^{(k)}\alpha(\mathbf{i}, \mathbf{j})} \tag{10}$$

i = 1, 2, …, n.

updated until it converges to the solution that optimizes an objective function. The objective

Algebraic reconstruction methods assume that the image is an array of unknowns, and the reconstruction problem can be established as a system of linear equations. The unknowns of this system are approximated with respect to the ray sums iteratively [17, 18]. In each iteration, current reconstruction is reprojected and updated according to how much it satisfies the

In the DT reconstruction, x represents the pixel values for x∈ ℜ<sup>N</sup>, whereas f represents the observed data for f ∈ ℜ<sup>M</sup>. The weighting value G ∈ ℜ<sup>M</sup>�<sup>N</sup> was created by using a forward projection, and the matrix G combines the submatrices of each projection. Gmn is defined as the length of intersection of the mth ray with the nth cell. DT reconstruction starts with an initial guess x<sup>0</sup> and computes x<sup>1</sup> by projecting x<sup>0</sup> onto the first hyperplane in Eq. 4. The algebraic

reconstruction technique (ART) update procedure (or error correction) is shown below:

þ β

measurement bi separately, that is, the kth iteration consists of a sweep through the m projection measurements. The algorithm iterates through the equations periodically; therefore,

The ART algorithm updates the image vector per ray such that the update satisfies only a single equation representing the corresponding ray integral. By contrast, the simultaneous IR technique (SIRT) algorithm updates the image vector after all equations are considered. The

The simultaneous ART (SART) algorithm [20] is proposed as a combination of the ART and SIRT algorithms. SART updates the superior implementation of ART and is based on a simultaneous update of the current reconstruction similar to SIRT. In the SART algorithm, the update procedure is applied for all rays in a given scan direction (projection) instead of each ray separately (similar to conventional ART) or instead of all rays simultaneously (similar to

Xm i aij bi � a ! i x !ð Þ <sup>k</sup>�<sup>1</sup>

P<sup>n</sup> <sup>j</sup>¼<sup>1</sup> aij

bi � a ! i x !ð Þ <sup>k</sup>�<sup>1</sup>

!¼ ð Þ ai1; ai2;…; ain is the <sup>i</sup>th row of the projection matrix, and <sup>β</sup> is the relaxation param-

a ! i � � � � � �

is the normalization factor. The update is performed for each projection

<sup>2</sup> a !

!ð Þ <sup>k</sup>�<sup>1</sup>

x !ð Þ<sup>k</sup> <sup>¼</sup> <sup>x</sup>

update procedure of SIRT is given in Eq. 6 according to [19].

!ð Þ <sup>k</sup>�<sup>1</sup>

<sup>þ</sup> <sup>β</sup> <sup>1</sup> P<sup>m</sup> <sup>i</sup> aij

x !ð Þ<sup>k</sup> <sup>¼</sup> <sup>x</sup>

SIRT). The SART update is expressed as follows:

G<sup>x</sup> ¼ f (4)

<sup>i</sup> (5)

(6)

function defines the criteria of the reconstruction solution.

projection measurements.

70 Medical Imaging and Image-Guided Interventions

where ai

eter (1.0). a

! i � � � � � � 2

i = (i-1) mod (m) + 1.

#### 2.4. Total variation minimization reconstruction algorithm

An iterative algorithm using total variation (TV)-based compressive sensing was recently developed for volume image reconstruction from a tomographic scan [21]. The image TV has been used as a penalty term in iterative image reconstruction algorithms [21]. The TV of an image is defined as the sum of the first-order derivative magnitudes for all pixels in the image. TV minimization is an image domain optimization method associated with compressed sensing theory [21]. As TV-minimization IR for image reconstruction, the TV-minimization IR technique is the SART [20] with algebraic IR for constraining the TV-minimization problem, which is called SART-TV [21]. In TV-minimization IR, adding a penalty to the data–fidelity–objective function tends to smooth out noise in the image while preserving edges within the image [21–25].

SART-TV minimizes the Rudin, Osher, and Fatemi (ROF) model [26], that is, SART-TV also takes into account the image information when minimizing the TV of the image. If only the TVminimization step was performed in the rest of the algorithms, the result would be a flat image; alternatively, the ROF model ensures that the image does not change considerably. The SART-TV optimal parameters include the iteration number for TV minimization [100 (np)] and the length of each gradient-descent step [50 (q)]. These optimal parameters have been shown to be very relevant to image quality [21, 24]. In our study, we minimized the image quality by using these optimal parameters for the SART-TV algorithms.

The SART-TV algorithm is expressed in the form of a pseudo code as follows:

$$
\beta = 1.0; \; \eta p; 100 \; \eta \text{-} 50
$$

$$
\overrightarrow{\text{x}} = 0
$$

repeat main loop

x0 !¼<sup>x</sup> !

(SART)

$$\text{for } i = 1: N\_{d\prime} \quad \overrightarrow{\mathbf{x}} = \overrightarrow{\mathbf{x}} + \beta \frac{1}{\sum\_{i \in \Omega\_{\boldsymbol{t}}} A\_{i\dot{\boldsymbol{\jmath}}}} \sum\_{\boldsymbol{i} \in \Omega\_{\boldsymbol{t}}} A\_{i\dot{\boldsymbol{\jmath}}} \frac{b\_{i} - \overrightarrow{A}\_{i} \cdot \overrightarrow{\mathbf{x}}}{\sum\_{j=1}^{n} A\_{i\dot{\boldsymbol{\jmath}}}} $$
 
$$\text{for } \dot{\boldsymbol{i}} = 1: N\_{\dot{\boldsymbol{\nu}}} \quad \text{if } \boldsymbol{\chi}\_{i} < 0 \text{ then } \boldsymbol{\chi}\_{i} = 0$$

(enforce positivity)

xres ! <sup>¼</sup> <sup>x</sup> !

4. Virtual monochromatic DT imaging using DE

combination of the attenuation coefficients of two basis materials [10]:

r 

1

<sup>μ</sup>ð Þ¼ <sup>r</sup>; <sup>E</sup> <sup>μ</sup>

Given that the photoelectric effect and Compton scattering of X-ray photons in the diagnostic range (E < 140 keV) are the predominant mechanisms responsible for X-ray attenuation (monochromatic X-ray), the mass attenuation coefficient of any material can be expressed with sufficient accuracy as a linear combination of the photoelectric and Compton attenuation coefficients. Consequently, the mass attenuation coefficient can also be expressed as a linear

ð Þ� <sup>E</sup> <sup>r</sup>1ð Þþ <sup>r</sup> <sup>μ</sup>

where the basis materials exhibit different photoelectric and Compton characteristics, (μ/r)i(E) and i = 1, 2 denote the mass attenuation coefficients of the two basis materials, and ri(r) and

In principle, DE can only accurately decompose a mixture into two materials. Therefore, for DE measurement–based mixture decomposition into three constitutive materials, a third constituent must be provided to solve for three unknowns by using only two spectral measurements. In one solution, the sum of the volumes of the three constituent materials is assumed to be

r 

2

) of the two basis materials at location r.

ð Þ� E r2ð Þr (11)

State-Of-The-Art X-Ray Digital Tomosynthesis Imaging http://dx.doi.org/10.5772/intechopen.81667 73

4.1. Theory of material decomposition processing

Figure 1. Illustration of the CD phantom used in this study.

i = 1, 2 denotes the local densities (g/cm<sup>3</sup>

for i = 1:np do TV-minimization loop

$$
\Delta \mathbf{x} = \left| \overrightarrow{\mathbf{x}} - \overrightarrow{\mathbf{x}}\_0 \right|
$$

$$
\overrightarrow{d} \mathbf{x} = \nabla\_{\overrightarrow{\mathbf{x}}} \left\| \overrightarrow{\mathbf{x}} \right\|\_{TV}
$$

$$
\hat{d} \mathbf{x} = \overrightarrow{d} \mathbf{x} / \left| \overrightarrow{d} \mathbf{x} \right|
$$

$$
\overrightarrow{\mathbf{x}} = \overrightarrow{\mathbf{x}} - q^\* \Delta \mathbf{x}^\* \hat{d} \mathbf{x}
$$

end

return x ! res
