Open access peer-reviewed chapter

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

By Tsutomu Gomi

Submitted: June 11th 2018Reviewed: September 25th 2018Published: November 5th 2018

DOI: 10.5772/intechopen.81667

Downloaded: 229

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.

Keywords

  • tomosynthesis
  • reconstruction
  • virtual monochromatic imaging
  • polychromatic imaging
  • dual-energy imaging

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 I0=δcdand that of the X-rays that passed through the structure at the location δcdas I=δcd. The image data dF=δcdwere calculated as follows:

dFδcd=lnI0δcdIδcdE1

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

fFDK=θθR2U2d˜FδcdE2
UXYδ=R+Xcosδ+Ysinδ
cXYδ=RXsinδ+YcosδR+Xcosδ+Ysinδ
dXYZδ=ZRR+Xcosδ+Ysinδ
d˜Fδcd=RR2+c2+d2dFδcdgPγ

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

gpγ=γsinα2E3

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.

Gx=fE4

In the DT reconstruction, xrepresents the pixel values for xN, whereas frepresents the observed data for fM. The weighting value GM×Nwas created by using a forward projection, and the matrix Gcombines the submatrices of each projection. Gmnis defined as the length of intersection of the mth ray with the nth cell. DT reconstruction starts with an initial guess x0 and computes x1 by projecting x0 onto the first hyperplane in Eq. 4. The algebraic reconstruction technique (ART) update procedure (or error correction) is shown below:

xk=xk1+βbiaixk1ai2aiE5

where ai=ai1ai2ainis the ith row of the projection matrix, and β is the relaxation parameter (1.0). ai2is 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].

xk=xk1+β1imaijimaijbiaixk1j=1naijE6

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:

xk=xk1+β1ΣiΩtaijiΩtaijbiaixk1j=1naijE7

where Ωtis 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 αijrepresents 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)]:

xk+1i=xkiCkiE8
Cki=1j=1nαijj=1Jyji=1nαijxkiαijE9
xk+1i=1j=1nαijxkij=1Jyjαiji=1nxkαijE10

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 (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:

β=1.0; np;100 q;50
x=0
repeat main loop
x0=x

(SART)

fori=1:Nd,x=x+β1iΩtAijiΩtAijbiAixj=1nAij
fori=1:Ni,ifxi<0thenxi=0

(enforce positivity)

xres = x
for i = 1:np do TV-minimization loop
Δx=xx0
dx=xxTV
d̂x=dx/dx
x=xqΔxd̂x
end
returnxres

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, CaCO3) 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]:

μrE=μρ1Eρ1r+μρ2Eρ2rE11

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/cm3) 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 (fn) of the CaCO3 (fcaco3, local density; 2.711 g/cm3), PMMA (fPMMA, local density; 1.17 g/cm3), and epoxy resin (fepoxy, local density; 1.11 g/cm3) in the phantom.

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

μrE=μρ1Eρ1r+μρ2Eρ2r+μρ3Eρ3rE12

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

IL=PLEexpμρ1EL1μρ2EL2μρ3EL3dEE13
IH=PHEexpμρ1EL1μρ2EL2μρ3EL3dEE14
L1L2L3=1.0E15
L1=ρ1rdl
L2=ρ2rdl
L3=ρ3rdl

where PL(E) represents the low-energy primary intensities, PH(E) represents the high-energy primary intensities, IL represents the low-energy attenuated intensities, and IH denotes the high-energy attenuated intensities. The equivalent densities (g/cm2; L1, L2, and L3) 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 L1, L2, and L3 are the unknown materials. Therefore, the basis material decomposition can be accomplished by solving simultaneous equations to calculate the values of L1, L2, and L3 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 CaCO3, epoxy resin, and PMMA into the XCOM program developed by Berger and Hubbell [27]. The curves show that the linear attenuation coefficient of CaCO3 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 CaCO3, epoxy resin image, and PMMA by using the following process.

Figure 2.

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

Eqs. (13) and (14) were used to calculate the values for IL_caco3, IL_PMMAr, IL_epoxy, IH_caco3, IH_PMMA, and IH_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):

fCaCO3fPMMAfepoxy=IL_CaCO3IL_PMMAIL_epoxyIH_CaCO3IH_PMMAIH_epoxy1.01.01.01DTELDTEH1.0E16
fCaCO3IL_CaCO3+fPMMAIL_PMMA+fepoxyIL_epoxy=DTEL
fCaCO3IH_CaCO3+fPMMAIH_PMMA+fepoxyIH_epoxy=DTEH
fCaCO3+fPMMA+fepoxy=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 CaCO3, epoxy, and PMMA.

Virtual monochromatic processing is performed according to Eq. 15:

Mono_p_img=fCaCO3μρCaCO3E+fPMMAμρPMMAE+fepoxyμρepoxyEE17

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 μm2). 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:

RMSE=i=1nykyk2nE18

where ykis the observed image [current reconstructed image (in-focus plane)], ykis 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:

SSIMxy=lxyαcxyβsxyγE19

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

α=β=γ=1.0

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

MSSIMXY=1Mj=1MSSIMxiyjE20

where X and Y are the reference [previous reconstructed image (in-focus plane)] and objective [current reconstructed image (in-focus plane)] images, respectively; xiand yjare 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 φ; CaCO3, 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:

CNR=μFeatureμBGσBGE21

where μFeatureis the mean object pixel value, μBGis the mean background area pixel value, and σBGis 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 CaCO3, 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 (CaCO3; 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 (CaCO3; 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.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tsutomu Gomi (November 5th 2018). State-Of-The-Art X-Ray Digital Tomosynthesis Imaging, Medical Imaging and Image-Guided Interventions, Ronnie A. Sebro, IntechOpen, DOI: 10.5772/intechopen.81667. Available from:

chapter statistics

229total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

A High-Performance System Architecture for Medical Imaging

By Tassadaq Hussain, Amna Haider, Muhammad Shafique and Abdelmalik Taleb-Ahmed

Related Book

First chapter

Longitudinal Changes of Structural and Functional Connectivity and Correlations with Neurocognitive Metrics

By Yongxia Zhou

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us