State-Of-The-Art X-Ray Digital Tomosynthesis Imaging

Abstract

Digital tomosynthesis (DT) is a notable modality in medical imaging because it shows the spread of the target area with lower radiation dose relative to computed tomography. In this section, we describe the technique in two parts: (1) image quality (contrast) and (2) DT image reconstruction algorithms, including state-of-the-art total variation minimization reconstruction algorithms with single-energy X-ray conventional polychromatic imaging and novel dual-energy (DE) virtual monochromatic imaging. The novel DE virtual monochromatic image-processing algorithm provides adequate overall performance (especially, reduction of beam-hardening, reduction of noise). The DE virtual monochromatic image-processing algorithm appears to be a promising new option for imaging in DT because it provides three-dimensional visualizations of high-contrast images that are far superior to those of images processed by using conventional single-energy polychromatic image-processing algorithms.

1. Introduction

Approximately 30 years have passed since Dr. Hounsfield developed a practical computed tomography (CT) system. The arrival of CT devices in the field of medical diagnosis has led to a revolution equivalent to the discovery of X-rays by Dr. Roentgen. Since then, most researchers who aim to improve medical diagnosis quality have worked toward the functional improvement of CT instruments, thus leading to the development of fan-beam CT, helical scan CT, multislice CT, and cone-beam CT instruments. These new instruments reduce the time needed for image reconstruction and significantly improve image quality. Given the increasing demand for better technology, there has been continued research and development of high-performance CT instruments.

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, 4, 5, 6, 7, 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 have been explored.

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, 11, 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.

2. DT reconstruction

Existing tomosynthesis algorithms can be divided into three categories: (1) backprojection (BP), (2) filtered backprojection (FBP), and (3) iterative reconstruction (IR) algorithms.

2.1. BP

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 parallel to the detector.

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_0=\left(\delta ,c,d\right)$ and that of the X-rays that passed through the structure at the location $\left(\delta ,c,d\right)$ as $I=\left(\delta ,c,d\right)$ . The image data $d^F=\left(\delta ,c,d\right)$ were calculated as follows:

$$d^F\left(\delta ,c,d\right)=\mathit{ln}\frac{I_0\left(\delta ,c,d\right)}{I\left(\delta ,c,d\right)}$$ E1

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

$$f^{\mathit{FDK}}={\int }_{-\theta }^{\theta }\frac{R^2}{U^2}{\tilde{d}}^F\left(\delta ,c,d\right)\mathit{D\delta }$$ E2

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

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

$$d\left(X,Y,Z,\delta \right)=Z\frac{R}{R+X\text{cos}\delta +\mathit{Ys}\text{in}\delta }$$

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

where the projection angle $\delta$ , fan-angle $\mathrm{\gamma }$ , source trajectory $R$ , acquisition angle $\theta$ , and $g^p\left(\gamma \right)$ represent a convolution with the Ramachandran–Lakshminarayanan filter. The Ramachandran–Lakshminarayanan (ramp) filter is shown below:

$$g^p\left(\gamma \right)={\left(\frac{\gamma }{\mathit{sin}}\right)}^2$$ E3

2.3. IR

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.

$$G_{\mathrm{x}}=f$$ E4

In the DT reconstruction, $\mathrm{x}$ represents the pixel values for $\mathrm{x}\in {\mathrm{\Re }}^N$ , whereas $f$ represents the observed data for $f\in {\mathrm{\Re }}^M$ . The weighting value $G\in {\mathrm{\Re }}^{M\times N}$ was created by using a forward projection, and the matrix $G$ combines the submatrices of each projection. $G_{\mathit{mn}}$ is defined as the length of intersection of the m^th ray with the n^th cell. DT reconstruction starts with an initial guess x⁰ and computes x¹ by projecting x⁰ onto the first hyperplane in Eq. 4. The algebraic reconstruction technique (ART) update procedure (or error correction) is shown below:

$${{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k\right)}={{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k-1\right)}+\beta \frac{b_i-{{\displaystyle \overrightarrow{a}}}_i{{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k-1\right)}}{{‖{{\displaystyle \overrightarrow{a}}}_i‖}^2}{{\displaystyle \overrightarrow{a}}}_i$$ E5

where ${\overrightarrow{a}}_i=\left(a_{i1},a_{i2},\dots ,a_{\mathit{\text{in}}}\right)$ is the ith row of the projection matrix, and β is the relaxation parameter (1.0). ${‖{{\displaystyle \overrightarrow{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].

$${{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k\right)}={{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k-1\right)}+\beta \frac{1}{{\sum }_i^m{a}_{\mathit{ij}}}\sum _i^m{a}_{\mathit{ij}}\frac{b_i-{{\displaystyle \overrightarrow{a}}}_i{{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k-1\right)}}{{\sum }_{j=1}^n{a}_{\mathit{ij}}}$$ E6

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:

$${{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k\right)}={{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k-1\right)}+\beta \frac{1}{{\mathrm{\Sigma }}_{i\in {\mathrm{\Omega }}_t}{a}_{\mathit{ij}}}\sum _{i\in {\mathrm{\Omega }}_t}{a}_{\mathit{ij}}\frac{b_i-{{\displaystyle \overrightarrow{a}}}_i{{\displaystyle \overrightarrow{\mathrm{x}}}}^{\left(k-1\right)}}{{\sum }_{j=1}^n{a}_{\mathit{ij}}}$$ E7

where ${\mathrm{\Omega }}_t$ 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 $\alpha \left(i,j\right)$ 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)]:

$${\mathrm{x}}^{\left(k+1\right)}\left(i\right)={\mathrm{x}}^{\left(k\right)}\left(i\right)C^{\left(k\right)}\left(i\right)$$ E8

$$C^{\left(k\right)}\left(i\right)=\frac{1}{{\sum }_{j=1}^n\alpha \left(i,j\right)}\sum _{j=1}^J\frac{y\left(j\right)}{{\sum }_{i=1}^n\alpha \left(i,j\right){\mathrm{x}}^{\left(k\right)}\left(i\right)}\alpha \left(i,j\right)$$ E9

$${\mathrm{x}}^{\left(k+1\right)}\left(i\right)=\frac{1}{{\sum }_{j=1}^n\alpha \left(i,j\right)}{\mathrm{x}}^{\left(k\right)}\left(i\right)\sum _{j=1}^J\frac{y\left(j\right)\alpha \left(i,j\right)}{{\sum }_{i=1}^n{\mathrm{x}}^{\left(k\right)}\alpha \left(i,j\right)}$$ E10

i = 1, 2, …, n.

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, 22, 23, 24, 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 TV-minimization 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 ( $\mathit{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; \mathit{np};100 q;50$$

$$\overrightarrow{x}=0$$

repeat main loop

$${\overrightarrow{x}}_0=\overrightarrow{x}$$

(SART)

$$\text{for}i=1:N_d,\overrightarrow{x}=\overrightarrow{x}+\beta \frac{1}{{\sum }_{i\in {\mathrm{\Omega }}_t}{A}_{\mathit{ij}}}\sum _{i\in {\mathrm{\Omega }}_t}{A}_{\mathit{ij}}\frac{b_i-{{\displaystyle \overrightarrow{A}}}_i\cdot {\displaystyle \overrightarrow{x}}}{{\sum }_{j=1}^n{A}_{\mathit{ij}}}$$

$$\text{for}i=1:N_i,\text{if}{x}_i<0\text{then}{x}_i=0$$

(enforce positivity)

$${\overrightarrow{x}}_{\mathit{res}} = \overrightarrow{x}$$

for i = 1:np do TV-minimization loop

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

$${\displaystyle \overrightarrow{d}}x={\nabla }_{{\displaystyle \overrightarrow{x}}}{‖{\displaystyle \overrightarrow{x}}‖}_{\mathit{TV}}$$

$${\displaystyle \widehat{d}}x={\displaystyle \overrightarrow{d}}x/\left|{\displaystyle \overrightarrow{d}}x\right|$$

$${\displaystyle \overrightarrow{x}}={\displaystyle \overrightarrow{x}}-q^{\ast }{\mathit{\Delta x}}^{\ast }{\displaystyle \widehat{d}}x$$

end

$$\text{return}{{\displaystyle \overrightarrow{x}}}_{\mathit{res}}$$

3. Phantom specification

For the evaluation of low-contrast resolution, a contrast-detail (CD) phantom was used with epoxy slabs. The CD phantoms of different diameters (signal region, CaCO₃) and thicknesses were arranged within the epoxy slabs. For X-ray imaging, we placed polymethyl methacrylate (PMMA) slabs (200 × 200 mm) with 50 mm thickness on the top of the CD phantom (Figure 1).

Figure 1.

Illustration of the CD phantom used in this study.

4. Virtual monochromatic DT imaging using DE

4.1. Theory of material decomposition processing

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 combination of the attenuation coefficients of two basis materials [10]:

$$\mu \left(r,E\right)={\left(\frac{\mu }{\rho }\right)}_1\left(E\right)\cdot {\rho }_1\left(r\right)+{\left(\frac{\mu }{\rho }\right)}_2\left(E\right)\cdot {\rho }_2\left(r\right)$$ E11

where the basis materials exhibit different photoelectric and Compton characteristics, (μ/ρ)_i(E) and i = 1, 2 denote the mass attenuation coefficients of the two basis materials, and ρ_i(r) and i = 1, 2 denotes the local densities (g/cm³) of the two basis materials at location r.

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 equivalent to the volume of the mixture (volume or mass conservation) [11]. In this work, we used a simple projection space (prereconstruction) decomposition method to estimate the material fractions (f_n) of the CaCO₃ (f_caco3, local density; 2.711 g/cm³), PMMA (f_PMMA, local density; 1.17 g/cm³), and epoxy resin (f_epoxy, local density; 1.11 g/cm³) in the phantom.

Three basis materials can also be expressed as a linear combination of the attenuation coefficients:

$$\mu \left(r,E\right)={\left(\frac{\mu }{\rho }\right)}_1\left(E\right)\cdot {\rho }_1\left(r\right)+{\left(\frac{\mu }{\rho }\right)}_2\left(E\right)\cdot {\rho }_2\left(r\right)+{\left(\frac{\mu }{\rho }\right)}_3\left(E\right)\cdot {\rho }_3\left(r\right)$$ E12

In DE acquisition, the detected image intensity can be expressed as follows:

$$I_L=\int P_L\left(E\right)\mathit{exp}\left\{-{\left(\frac{\mu }{\rho }\right)}_1\left(E\right)\cdot L_1-{\left(\frac{\mu }{\rho }\right)}_2\left(E\right)\cdot L_2-{\left(\frac{\mu }{\rho }\right)}_3\left(E\right)\cdot L_3\right\}\mathit{dE}$$ E13

$$I_H=\int P_H\left(E\right)\mathit{exp}\left\{-{\left(\frac{\mu }{\rho }\right)}_1\left(E\right)\cdot L_1-{\left(\frac{\mu }{\rho }\right)}_2\left(E\right)\cdot L_2-{\left(\frac{\mu }{\rho }\right)}_3\left(E\right)\cdot L_3\right\}\mathit{dE}$$ E14

$${L_1}^{\ast }{L_2}^{\ast }{L}_3=1.0$$ E15

$$L_1=\int {\rho }_1\left(r\right)\mathit{dl}$$

$$L_2=\int {\rho }_2\left(r\right)\mathit{dl}$$

$$L_3=\int {\rho }_3\left(r\right)\mathit{dl}$$

where P_L(E) represents the low-energy primary intensities, P_H(E) represents the high-energy primary intensities, I_L represents the low-energy attenuated intensities, and I_H denotes the high-energy attenuated intensities. The equivalent densities (g/cm²; L₁, L₂, and L₃) of the three basis materials must be determined for each ray path. Eqs. (13, 14, 15) can be solved for the equivalent area density, where L₁, L₂, and L₃ are the unknown materials. Therefore, the basis material decomposition can be accomplished by solving simultaneous equations to calculate the values of L₁, L₂, and L₃ from the measured projection pixel values [12]. By using the density corresponding to the area with the three basis materials, the linear attenuation coefficient μ(r, E) can be calculated for any photon.

We used the local and area densities for each material to calculate the theoretical linear attenuation coefficient curves shown in Figure 2; these values were generated by inputting the chemical compositions of the CaCO₃, epoxy resin, and PMMA into the XCOM program developed by Berger and Hubbell [27]. The curves show that the linear attenuation coefficient of CaCO₃ decreases more rapidly than those of the foam epoxy and PMMA in the energy band <100 keV. Finally, we used the projection space decomposition approach to generate material decomposition images for the CaCO₃, epoxy resin image, and PMMA by using the following process.

Figure 2.

The linear attenuation coefficients of CaCO₃, epoxy, and PMMA with respect to the photons.

Eqs. (13) and (14) were used to calculate the values for I_L__caco3, I_L__PMMAr, I_L__epoxy, I_H__caco3, I_H__PMMA, and I_H__epoxy as simulated attenuation intensities of these materials at the two energy levels. These values were then used to construct a sensitivity matrix, and the material fractions were obtained from the inverse of this matrix (Eq. 16):

$$\left[\begin{array}{l}{f}_{\mathit{\text{CaCO}}3}\\ f_{\mathit{\text{PMMA}}}\\ f_{\mathit{\text{epoxy}}}\end{array}\right]={\left[\begin{array}{lll}{I}_{L\_\mathit{\text{CaCO}}3}& I_{L\_\mathit{\text{PMMA}}}& I_{L\_\mathit{\text{epoxy}}}\\ I_{H\_\mathit{\text{CaCO}}3}& I_{H\_\mathit{\text{PMMA}}}& I_{H\_\mathit{\text{epoxy}}}\\ 1.0& 1.0& 1.0\end{array}\right]}^{-1}\left[\begin{array}{l}{\mathit{DT}}_{\mathit{EL}}\\ {\mathit{DT}}_{\mathit{EH}}\\ 1.0\end{array}\right]$$ E16

$$f_{\mathit{\text{CaCO}}3}{I}_{L\_\mathit{\text{CaCO}}3}+f_{\mathit{\text{PMMA}}}{I}_{L\_\mathit{\text{PMMA}}}+f_{\mathit{\text{epoxy}}}{I}_{L\_\mathit{\text{epoxy}}}={\mathit{DT}}_{\mathit{EL}}$$

$$f_{\mathit{\text{CaCO}}3}{I}_{H\_\mathit{\text{CaCO}}3}+f_{\mathit{\text{PMMA}}}{I}_{H\_\mathit{\text{PMMA}}}+f_{\mathit{\text{epoxy}}}{I}_{H\_\mathit{\text{epoxy}}}={\mathit{DT}}_{\mathit{EH}}$$

$$f_{\mathit{\text{CaCO}}3}+f_{\mathit{\text{PMMA}}}+f_{\mathit{\text{epoxy}}}=1.0$$

4.2. Virtual monochromatic image processing

After decomposition by matrix inversion, the “inv” function available in MATLAB was used (Mathworks; Natick, MA, USA); this function constrains the possible fraction to [0,1] while imposing a sum of one. Accordingly, the processing pipeline yields three material fraction outputs corresponding to CaCO₃, epoxy, and PMMA.

Virtual monochromatic processing is performed according to Eq. 15:

$$\mathit{\text{Mono}}\_p\_\mathit{img}={f_{\mathit{\text{CaCO}}3}}^{\ast }{\left(\frac{\mu }{\rho }\right)}_{\mathit{\text{CaCO}}3}\left(E\right)+{f_{\mathit{\text{PMMA}}}}^{\ast }{\left(\frac{\mu }{\rho }\right)}_{\mathit{\text{PMMA}}}\left(E\right)+{f_{\mathit{\text{epoxy}}}}^{\ast }{\left(\frac{\mu }{\rho }\right)}_{\mathit{\text{epoxy}}}\left(E\right)$$ E17

where Mono_p_img is the virtual monochromatic projection image, and [μ/ρ]_caco3(E), [μ/ρ]_PMMA(E), and [μ/ρ]_epoxy(E) are the mass attenuation coefficients of each material. The generated virtual monochromatic X-ray projection image was reconstructed by using each algorithm for energies of 60, 80, 100, 120, and 140 keV. The real projection data acquired on a DT system were used for reconstruction. All image reconstruction calculations, including DE material decomposition processing and reconstruction, as well as FBP, SART, SART-TV, virtual monochromatic processing, and MLEM, were implemented in MATLAB.

5. DT system

5.1. DT overview

The DT system (SonialVision Safire II; Shimadzu Co., Kyoto, Japan) comprised an X-ray tube [anode: tungsten with rhenium and molybdenum; real filter: inherent; aluminum (1.1 mm), additional; aluminum (0.9 mm), and copper (0.1 mm)] with a 0.4 mm focal spot and 362.88 × 362.88 mm amorphous selenium digital flat-panel detector (detector element, 150 × 150 μm²). The source-to-isocenter and isocenter-to-detector distances were 924 and 1100 mm, respectively (antiscatter grid, focused type; grid ratio, 12:1).

5.2. DE-DT acquisition

Collimator motion was synchronized by measuring the misalignment of the low-voltage (60 kV) and high-voltage (120 kV) images at a constant tube motion. A large energy gap between low and high tube potential kVp imaging yields better material decomposition [28, 29, 30, 31, 32, 33]. We selected the abovementioned kV values because this study aimed to improve contrast and artifact reduction during DT acquisition while maintaining an imaging performance similar to that of conventional DT.

Pulsed X-ray exposures and rapid switching between low and high tube potential kVp values were used for DE-DT imaging. Linear system movement and a swing angle of 40° were used when performing tomography, and 37 low- and high-voltage projection images were sampled during a single tomographic pass. We used a low voltage, and each projection image was acquired at 416 mA. A 9.4 ms exposure time was used for low-voltage (60 kV) X-rays at 416 mA, and a 2.5 ms exposure time was used for high-voltage (120 kV) X-rays. To generate reconstructed tomograms of the desired height, we used a 768 × 7684 matrix with 32 bits (single-precision floating number) per image (pixel size, 0.252 mm/pixel; reconstruction interval, 1 mm).

6. Optimization parameter

The experiments were performed according to the scheme shown in Figure 3. A range of optional parameters have been identified for IR algorithms. Among these parameters, some are important for determining algorithmic behavior. In this study, we compared the root-mean-square error (RMSE) and mean structural similarity (MSSIM; reconstructed volume image from the previous iterations between the current iteration) to optimize the iteration numbers (i).

Figure 3.

For DT acquisition, the phantom was arranged parallel to the x–y detector plane.

The RMSE was defined in this study as follows:

$$\mathit{\text{RMSE}}=\sqrt{\frac{{\sum }_{i=1}^n{\left({{\displaystyle \stackrel{\wedge }{y}}}_k-y_k\right)}^2}{n}}$$ E18

where $y_k$ is the observed image [current reconstructed image (in-focus plane)], ${{\displaystyle \stackrel{\wedge }{y}}}_k$ is the referenced image [previous reconstructed image (in-focus plane)], and n is the number of compounds in the analyzed set.

The MSSIM of local patterns of luminance- and contrast-normalized pixel intensity were compared to determine the structural similarity (SSIM) index of contrast preservation. This image quality metric is based on the assumed suitability of the human visual system for extracting structure-based information [34].

The SSIM index between pixel values x and y was calculated as follows:

$$\mathit{\text{SSIM}}\left(x,y\right)={\left[l\left(x,y\right)\right]}^{\alpha }\cdot {\left[c\left(x,y\right)\right]}^{\beta }\cdot {\left[s\left(x,y\right)\right]}^{\gamma }$$ E19

where l is the luminance, c is the contrast, and s is the structure. Subsequently,

$$\alpha =\beta =\gamma =1.0$$

The MSSIM was then used to evaluate the overall image quality:

$$\mathit{\text{MSSIM}}\left(X,Y\right)=\frac{1}{M}\sum _{j=1}^M\mathit{\text{SSIM}}\left(x_i,y_j\right)$$ E20

where X and Y are the reference [previous reconstructed image (in-focus plane)] and objective [current reconstructed image (in-focus plane)] images, respectively; $x_i$ and $y_j$ are the image contents at the jth pixel; and M is the number of pixels in the image.

It was feasible to maintain a steady convergence of polychromatic IR images for inconsistency after the fourth iteration and a convergence of monochromatic IR images for inconsistency after the fifth iteration (Figures 4 and 5).

Figure 4.

The RMSE characteristics caused by the differences in the number of iterations of each IR algorithm. (60 and 120kVp; polychromatic, 60, 80, 100, 120, and 140 keV; virtual monochromatic).

Figure 5.

The MSSIM characteristics caused by differences in the number of iterations of each IR algorithm. (60 and 120kVp; polychromatic, 60, 80, 100, 120, and 140 keV; virtual monochromatic).

7. Evaluation

The contrast derived from the contrast-to-noise ratio (CNR) in the in-focus plane (15 mm φ; CaCO₃, 175 mg/ml, and 100 mg/ml) was also evaluated as a quantitative measure of the reconstructed image quality. In DT, the CNR is frequently used to estimate low-contrast detectability and was defined in this study as follows:

$$\mathit{CNR}=\frac{{\mu }_{\mathit{\text{Feature}}}-{\mu }_{\mathit{BG}}}{{\sigma }_{\mathit{BG}}}$$ E21

where ${\mu }_{\mathit{\text{Feature}}}$ is the mean object pixel value, ${\mu }_{\mathit{BG}}$ is the mean background area pixel value, and ${\sigma }_{\mathit{BG}}$ is the standard deviation of the background pixel values. The latter parameter includes the photon statistics, the electronic noise from the results, and the structural noise that could obscure the object. The sizes of all regions of interest (ROIs) used to measure the CNR were adjusted to an internal signal (ROI diameter; eight pixels).

8. Results

Figure 6 shows each density projection image generated by the material decomposition process. A novel DE virtual monochromatic image processing is performed from a density projection image, and Figures 7 and 8 show the image generated by each reconstruction algorithm. In the novel DE virtual monochromatic image processing performed according to Eq. 17, an image with decreased noise can be obtained as the photon energy was increased. In the conventional polychromatic image, noise is reduced in the high-voltage image. With regard to the contrast extraction capability of high-density objects, the novel DE virtual monochromatic images showed high CNRs. In particular, the SART-TV algorithm provided high-contrast extraction capability (Figure 9). Along with the decrease in the concentration of the object, the contrast extraction ability was almost the same for both monochromatic and virtual monochromatic imaging (Figure 10).

Figure 6.

DE material decomposition straight projection images. The DE material decomposition projections for CaCO₃, epoxy, and PMMA were window level = 0.59 and window width = 0.26, window level = 0.69 and window width = 0.33, and window level = 0.21 and window width = 0.24, respectively.

Figure 7.

Comparisons among the polychromatic (P) and virtual monochromatic (VM) images acquired by using the FBP reconstruction algorithm (window level = 0.05, window width = 0.11).

Figure 8.

Comparisons among the polychromatic (P) and virtual monochromatic (VM) images with each IR algorithm (P MLEM: Window level = 0.32, window width = 0.17; VM MLEM: Window level = 0.25, window width = 0.12; P&VM SART: Window level = 0.22, window width = 0.37; P&VM SART-TV: Window level = 0.24, window width = 0.34).

Figure 9.

Comparisons of the CNR of the in-focus plane images obtained by using each reconstruction algorithm with polychromatic and virtual monochromatic image processing (CaCO₃; 15 mm ϕ, 175 mg/ml).

Figure 10.

Comparisons of the CNR of in-focus plane images obtained by using each reconstruction algorithm with polychromatic and virtual monochromatic image processing (CaCO₃; 15 mm ϕ, 100 mg/ml).

9. Conclusions

This chapter showed that the novel DE virtual monochromatic image-processing algorithm yielded adequate overall performance. The novel DE virtual monochromatic images produced by using this algorithm yielded good results independent of the type of contrast present in the CD phantom. Furthermore, this processing algorithm successfully removed noise from the images, particularly at high contrast with high-density objects.

In summary, this DE virtual monochromatic image-processing algorithm appears to be a promising new option for DT imaging, as evidenced by the 3D visualizations of high-contrast images that were far superior to those of images processed by conventional single-energy polychromatic image-processing algorithms. The flexibility in the choice of imaging parameters, which is based on the desired final images and DT imaging conditions, of this novel DE virtual monochromatic image-processing algorithm may be beneficial to users.

Acknowledgments

We wish to thank Mr. Kazuaki Suwa and Yuuki Watanabe at Department of Radiology Dokkyo Medical University Koshigaya Hospital for support on experiment.

Conflict of interest

The authors declare that they have no conflict of interest.

Ethical approval

This article does not contain any studies with human participants or animals performed by any of the authors.

Informed consent

This articles does not contain patient data.