|Home | About | Journals | Submit | Contact Us | Français|
Digital tomosynthesis mammography (DTM) is a promising new modality for breast cancer detection. In DTM, projection-view images are acquired at a limited number of angles over a limited angular range and the imaged volume is reconstructed from the two-dimensional projections, thus providing three-dimensional structural information of the breast tissue. In this work, we investigated three representative reconstruction methods for this limited-angle cone-beam tomographic problem, including the backprojection (BP) method, the simultaneous algebraic reconstruction technique (SART) and the maximum likelihood method with the convex algorithm (ML-convex). The SART and ML-convex methods were both initialized with BP results to achieve efficient reconstruction. A second generation GE prototype tomosynthesis mammography system with a stationary digital detector was used for image acquisition. Projection-view images were acquired from 21 angles in 3° increments over a ±30° angular range. We used an American College of Radiology phantom and designed three additional phantoms to evaluate the image quality and reconstruction artifacts. In addition to visual comparison of the reconstructed images of different phantom sets, we employed the contrast-to-noise ratio (CNR), a line profile of features, an artifact spread function (ASF), a relative noise power spectrum (NPS), and a line object spread function (LOSF) to quantitatively evaluate the reconstruction results. It was found that for the phantoms with homogeneous background, the BP method resulted in less noisy tomosynthesized images and higher CNR values for masses than the SART and ML-convex methods. However, the two iterative methods provided greater contrast enhancement for both masses and calcification, sharper LOSF, and reduced inter-plane blurring and artifacts with better ASF behaviors for masses. For a contrast-detail phantom with heterogeneous tissue-mimicking background, the BP method had strong blurring artifacts along the x-ray source motion direction that obscured the contrast-detail objects, while the other two methods can remove the superimposed breast structures and significantly improve object conspicuity. With a properly selected relaxation parameter, the SART method with one iteration can provide tomosynthesized images comparable to those obtained from the ML-convex method with seven iterations, when BP results were used as initialization for both methods.
Mammography is the only proven method currently used for breast cancer screening. In mammography, the x rays transmitted through the compressed breast are recorded on a screen-film system (conventional mammography) or by a digital detector (digital mammography). For screening mammography, x-ray projection images of the breast are taken from two views, namely, the cranial-caudal (CC) view and the mediolateral-oblique (MLO) view. Other special views may be obtained during diagnostic work-up.
The CC-view and MLO-view images only contain two-dimensional (2D) projection information of three-dimensional (3D) anatomical structures, and thus the accuracy of breast cancer detection based on screening mammography is affected by both the fact that the breast cancer is often obscured by overlapping breast tissues, which may cause false-negative diagnoses, and the fact that the superimposed normal tissues mimic masses, which may cause false-positive diagnoses.
Digital tomosynthesis mammography (DTM) is a promising technique that can provide 3D structural information by reconstructing the whole imaged volume from a sequence of projection-view mammograms.1–4 It has been demonstrated.5,6 that DTM can reduce the camouflaging effect of the overlapping fibroglandular breast tissue, thus improving the conspicuity of subtle lesions. Several manufacturers of digital mammography systems have developed prototype DTM systems and are conducting pilot clinical trials to evaluate the utility of DTM.7,8
The concept of conventional tomosynthesis was introduced by Ziedses des Plantes.9 who published the first work of geometric tomography. Early works of three-dimensional geometric tomography system include those by Garrison et al.,10 Richards et al.,11 Miller et al.,12 and Grant,13 in which Grant introduced the term “tomosynthesis.” In early tomosynthesis systems, a discrete series of projection radio-graphs were acquired by moving the tube and film simultaneously in opposite directions. The reconstruction process involved simply shifting and adding the projection images taken at different projection angles. If the amount of shift between the projection images is chosen appropriately, a certain plane inside the imaged subject will be in focus. All features on this focused plane will be enhanced while those on other planes blurred. Recent advent of digital x-ray detectors provides the possibility of fast acquisition of multiple low-dose images of the subject.14,15 and digital reconstruction of the tomosynthesized slices. Thus, improved reconstruction quality can be achieved by using advanced digital reconstruction algorithms. In breast tomosynthesis reconstruction, the projection images are generally acquired at a limited number of views in a limited angular range. The number of images acquired is limited by the total dose which should be comparable to that used in conventional mammography. The reconstruction of the 3D volume of the compressed breast from 2D projection images represents a limited-angle cone-beam tomographic problem. Existing reconstruction methods for such a problem can be coarsely divided into four categories: backprojection algorithms, transform algorithms, algebraic reconstruction techniques, and statistical reconstruction algorithms. A review of tomosynthesis reconstruction techniques can be found in Dobbins et al.16
To explore the feasibility of the limited-angle cone-beam reconstruction methods for breast tomosynthesis, we selected several representative algorithms and applied them to sets of breast phantom data. The selected algorithms in this study included the back-projection (BP) method, the simultaneous algebraic reconstruction technique (SART), and the maximum likelihood method with the convex algorithm (ML-convex). We will discuss these methods in Sec. II.
A second generation GE prototype tomosynthesis mammography system in the breast imaging research laboratory at the University of Michigan Hospital was used in this study. Projection-view images were acquired from 21 angles by automatically moving the x-ray tube in 3° increments over a ±30° angular range in less than 8 s. In contrast, the first generation GE prototype tomosynthesis system, developed by GE in collaboration with the Breast Imaging Group at Massachusetts General Hospital,2 acquired projection-view images from 11 angles over a ±25° angular range in 5° increments in 7 s.
We used an American College of Radiology (ACR) phantom and designed three additional phantoms in this study. We investigated the dependence of the reconstructed image quality and artifacts of the SART and ML-convex methods on the number of iterations. To quantitatively evaluate the contrast, sharpness, and artifacts we employed the contrast-to-noise ratio (CNR), a line profile of features, and an artifact spread function (ASF). The noise behavior and the relative sharpness of these reconstruction algorithms were compared by using the noise power spectrum (NPS) and line profiles of thin wires.
As described in the Sec. I, existing reconstruction methods for digital tomosynthesis, or equivalently limited-angle cone-beam x-ray tomography, can be classified into four categories, including backprojection algorithms, transform algorithms, algebraic reconstruction techniques, and statistical reconstruction algorithms. We will discuss the principle of methods in each category in this section and introduce the methods investigated in this study.
Let the whole image volume be subdivided into J voxels, and the linear attenuation coefficient for the jth voxel be denoted by xj, 1 ≤ j ≤ J; the digital area detector contains I elements, and the ith ray, 1 ≤ i ≤ I, is defined as a line segment starting from the point x-ray source location to the center of the ith detector element. The number of rays is equal to the number of detector elements assuming one ray is traced for each element. The path length of the ith ray going through the jth voxel in the nth x-ray tube location (projection view) is denoted by aij,n, resulting in a matrix-vector form of the projection model as
where An is the projection matrix for the nth projection view with aij,n as its (i, j)th element; yn is the corresponding vector of the projection data, which can be derived from the pixel values of the detected image, 1 ≤ n ≤ N, where N is the total number of projection views. The ith projection value, yi,n, is proportional to the logarithmic transform of the ratio of the incident intensity (Io,n) and the transmitted intensity (Ii,n) of the ith ray
We stack the projection model in Eq. (1) for all projection views together as
This is the linear system model for the breast tomosynthesis reconstruction, and it is the basis for the backprojection method and algebraic reconstruction techniques.
In the BP method, the linear attenuation coefficient of a given voxel is estimated by averaging all pixel values corresponding to x rays traveling through this voxel. Specifically, assuming that there are a total of M rays going through the jth voxel over all projection views, the linear attenuation coefficient is estimated by
where the normalization factor lm is the total path-length of the mth ray intersecting with the whole imaged volume, i.e., the sum of the corresponding row in the projection matrix.
We choose the backprojection method in this comparative study as it is easy to implement and its result can be used as initial values for more sophisticated iterative methods to achieve efficient reconstruction. The shift-and-add method used in the original tomosynthesis.1,13 (or referred to as the multiple projection method)17 is a simplified version of the BP method. Some additional strategies have been employed with the BP method, including the order statistics operator18,19 and reproduction of artifact image,20 to reduce the interplane artifacts, which are typically very strong in BP reconstructed images, especially those from high-contrast features.
Transform algorithms include the filtered backprojection (FBP) method and other transfer function methods. In the FBP method,21,22 the Fourier slice theorem plays a fundamental role and projection images are transformed to the spatial frequency domain. With a parallel-beam approximation, the 2D Fourier transform of each projection image corresponds to one slice sample in the spatial frequency domain at this angle and thus multiple projection images can be combined to obtain a discrete sampling of the whole imaged volume. Wedge filters, combining with other low-pass filters, are used in the spatial frequency domain to weight each projection in order to compensate for incomplete and/or non-uniform sampling in the spatial frequency domain and to suppress high frequency components. 2D FBP has been extensively used in computed tomography (CT) reconstruction in which a large number of projection images are acquired over a full angular range, so the corresponding spatial frequency domain is well sampled. For 3D cone-beam breast tomosynthesis, the Feldkamp or FDK algorithm,23 an approximation of the cone-beam FBP algorithm, had been investigated but the reconstruction results were found to be very noisy and the structure details were not visible.24 Nevertheless it provided good-quality reconstruction in volumetric CT breast imaging.25 A modified FBP method has been used in a prototype breast DTM system for tomosynthesis reconstruction but no details of the method has been reported.26
In other transfer function methods, impulse response functions or point spread functions are specially designed based on the imaging geometry and the application. The convolution process again becomes a simple multiplication as the Fourier transform is applied to each imaged slice or the whole imaged volume. Inversion filtering is employed to restore the whole imaged volume and to suppress the inter-plane blurring either in a continuous function form;27,28 or in a matrix-vector form if discrete voxel arrays are reconstructed.16,29 Promising results have been obtained for breast tomosynthesis with the former approach using some specially designed inverse filters.28,30 The latter is also referred to as matrix inversion tomosynthesis (MITS) and has been successfully applied to hand joint and chest radiographs.16,29,31 We did not choose the FBP algorithm in this study since FBP is restricted by the imaging geometry, and for a very limited-angle 3D reconstruction problem as in breast tomosynthesis where the spatial frequency domain is far from completely sampled, the image quality of the reconstructed DTMs depends strongly on the filter design.
In algebraic reconstruction techniques,21,32,33 the tomographic inverse problem is formulated as solving a large-scale system of linear equations, as shown in Eq. (3), where each element of the coefficient matrix A contains the intersection path length of one x ray within one specific voxel. The reconstruction is accomplished by iteratively updating the unknown linear attenuation coefficients by minimizing the error between the measured and the calculated projection data. For clarity and without loss of generality, we define a complete iteration of algebraic reconstruction methods as when all projection data over all projection views have been processed exactly once.
The original method in this family of algebraic reconstruction techniques is known by the same name—ART. In ART, the linear attenuation coefficients are updated in a “ray-by- ray” manner, i.e., updated after each “ray” equation (one row) of Eq. (3) is processed. All voxels along the ray under consideration are updated by the difference between the detected and computed pixel value. This difference is backprojected along the ray and contributes to each voxel in proportion to the path length of the ray inside this voxel.
Since only a single projection value is used to update the linear attenuation coefficients at a time, ART has fast convergence speed but will converge to a least squares solution which can be very noisy for severely ill-posed inverse problems such as limited-angle tomosynthesis reconstruction. To improve the ART method, variations on its implementation have been proposed. Depending on the different amount of projection data and the method used to update the current estimation, ART has been modified to other methods such as the SART34,35 and the simultaneous iterative reconstruction technique (SIRT).21,32 In SART, the linear attenuation coefficient of each voxel will not be updated until after all rays in one projection view have been processed once; while in SIRT, the update is performed after all projection views have been processed. Algebraic reconstruction techniques are flexible in regard to the measurement geometry and thus well suited for tomosynthesis reconstruction. Some studies have been conducted on applying the algebraic reconstruction techniques to tomosynthesis problem,36,37 but no application has been published on breast imaging to date.
In this study, we chose the SART method because this block-action strategy is a good trade-off between ART and SIRT. In ART, the measurement noise could be strongly amplified due to the frequent updating of the unknowns and the intrinsic ill-posed attribute of the inverse problem; whereas in SIRT, the convergence speed will be slow because the update is averaged over all rays in all projection views and the reconstruction could be overly smoothed. Breast tomosynthesis deals with a very large-scale linear system. To reconstruct efficiently, the chosen algorithm should converge to an acceptable solution within a small number of iterations.
As discussed above, for SART method, an update is performed after all rays in one projection view have been processed. The number of updates in one complete iteration of SART is equal to the number of projection views N. Let the linear attenuation coefficient in the jth voxel of the imaged volume at the completion of the kth iteration be denoted as . For the (k+1)th iteration, we first define the initial linear attenuation coefficient values as . Assume that there are a total of Mn rays in the nth projection view going through the jth voxel. The index of this subset of the I rays in the nth projection view is denoted by m, where 1 ≤ m ≤ Mn. The update of the linear attenuation coefficient at the jth voxel during processing of the nth projection view is defined as
for 1 ≤ j ≤ J. After all N projection views have been processed, the resulting estimate is given by for 1 ≤ j ≤ J.
For iterative algebraic reconstruction techniques like SART, several important aspects need to be considered for practical implementation. First, the initial values for the iterative process should be chosen properly. Different choices have been suggested, including constant distribution of zero (or very small positive) values, constant distribution of averaged attenuation coefficient,33 or reconstruction results from other methods such as the backprojection method.33 In this study, we chose the voxel values from the BP reconstruction as initial distributions for both the SART and ML methods. Second, one needs to select an appropriate step size or relaxation parameter λ. Although some previous works have suggested using eigenvalues of the forward matrix A and the expected number of iterations together to set the λ value,33 it is not feasible to estimate the eigenvalues for a very large-scale linear system such as that used in DTM. We therefore used a simple strategy in which λ is set to decrease over a very limited number of iterations. Third, it has been reported that the access order of projection views has a strong effect on the practical performance of algebraic reconstruction techniques in case of a full, or almost full, range of projections,38–40 except for SIRT method which uses the whole set of projection views once in each iteration. It has been suggested that the order of processed projections should be arranged such that successive projections are well “separated,” i.e., having least correlation in term of information. For DTM, generally only a limited number of projection views are acquired in a limited angular range, changing the access order of the projection views thus may not have strong effects on the reconstruction quality. We use a sequential access order in this study and leave the investigation of different access strategies to future studies. Last but not least, an appropriate stopping criterion of the iterative process has to be chosen. Previous studies have shown that early stopping can regularize ill-posed inverse problems. In practice, it is difficult to establish a specific criterion to stop the iteration automatically for a given tomosynthesis case. The number of iterations and the relaxation parameter are often predetermined based on both quantitative analysis of image quality measures, such as those investigated in this study, and visual comparison of reconstructed image quality for a given type of cases. The chosen parameters are then fixed and applied to this type of cases in future reconstructions.
All reconstruction algorithms discussed thus far are deterministic methods. On the other hand, statistical reconstruction methods tackle the inverse problem from a statistical point of view, assuming the unknown xj as a random variable following some specific probability distribution functions. Statistical reconstruction methods explicitly take into account the measurement statistics and noise model. In addition, similar to algebraic reconstruction techniques, statistical methods are suitable for any geometrical model.
Assume that the incident and transmitted x-ray intensities follow Poisson statistics and the measured intensity at different detector elements are statistically independent, the likelihood function or the conditional probability distribution of the measured intensities given the distribution of linear attenuation coefficients in the imaged volume, P(I|x), will have the computed detected intensity as the expectation value, as
where Ī i,n is the computed detected intensity at the ith detector element for the nth projection view
and is the total attenuation along the ith ray in the nth projection view or the ith component in the vector Anx. The likelihood function is defined over all projection views. The ML method is then used to estimate the unknown linear attenuation coefficients by searching for values such that the measured data has the largest probability to be obtained. When the components independent of the unknowns are ignored, the logarithm likelihood function can be written as
Due to the fact that there is no analytic solution for Eq. (8) in transmission tomography,41 it is difficult to search the entire space of unknown sets to find the ML solution. Various algorithms have been investigated to maximize the log-likelihood function iteratively, such as the expectation-maximization (EM) algorithm, convex algorithm, and gradient algorithm.42 Among them, the EM algorithm for transmission tomography involves calculation of the exponential components many times and thus is not efficient to implement. It has also been reported that the EM algorithm has a very slow convergence speed. In this work, we employed the convex algorithm (referred to as ML-convex method in this paper) to update the unknowns. The ML-convex method has been applied to breast tomosynthesis reconstruction by Wu et al.2,24
Using the ML-convex method, it can be derived that the update is given by
for 1 ≤ j ≤ J. Precisely, for each voxel, the update is proportional to the difference between the measured and computed intensities of the ith ray, weighted by the path length and normalized by the corresponding factors. More details about the convex algorithm for the ML method can be found in the work of Lange and Fessler42 and references therein.
Similar to the algebraic reconstruction techniques, the initial values and the step size are two important factors for the ML-convex algorithm, especially in practical implementation where an acceptable reconstruction result is expected to be accomplished within a small number of iterations. In this work, we used BP results as initial values, the same as those we used in the SART method, and kept the step size a constant value.
The imaging geometry of the second generation GE prototype digital tomosynthesis system for breast imaging research is shown in Fig. 1. The system has a CsI phosphor/a:Si active matrix flat panel digital detector with a pixel size of 0.1 mm × 0.1 mm and the raw image data are 16 bits. For tomosynthesis imaging, the x-ray tube is automatically rotated in 3° increments to acquire projection images at 21 different angles over a ±30° angular range in less than 8 s. The digital detector is stationary during image acquisition. The DTM system uses an Rh-target/Rh-filter x-ray source for all breast thicknesses. The kilovoltage setting ranges from 26 to 33 kVp and the total mA s of the 21 projection views ranges from 44 to 150 mA s. The imaging technique used depends on the breast thickness. No antiscatter grid is used. A single-projection image can be taken separately with the x-ray source fixed at 0° when the system is operated as a regular digital mammography system. For the ACR phantom, this system acquires a DTM using an exposure technique of Rh/Rh at 29 kV, the mean glandular dose of which is estimated to be about 250 mrad. Since the mean glandular dose limit for an ACR phantom imaged with a conventional mammography system is 200 mrad in the State of Michigan, the total glandular dose for the DTM is about 1.5 times of that for a conventional film mammogram. The output projection-view images are corrected for detector artifacts before being used for image reconstruction.
As shown in Fig. 1, the distance from the x-ray focal spot to the center of rotation is 64 cm and the plane along which the x-ray source rotates is perpendicular to the detector surface at the chest wall edge. The focal-spot-to-detector distance is 66 cm. The predefined imaged volume is 14.00 cm × 23.04 cm in area and 5 cm in thickness. The imaged volume thickness was set to be the maximum phantom thickness of 5 cm used in this study. Note that the actual detector size is 19.20 cm × 23.04 cm. We trim the output images to a smaller size that is slightly larger than the phantom dimensions to reduce reconstruction time. The imaged volume contains the whole breast phantoms. We subdivided the imaged volume into a set of voxels, of which the X and Y dimensions were set to be the same as the pixel size of the detector (0.1 mm × 0.1 mm) while the Z dimension (the slice thickness) was set to 1 mm.
For the projection model, we have developed an algorithm to calculate the path length of a primary x-ray intersecting each voxel within the breast volume. For a given comparison, all reconstruction algorithms used the same set of projection view images as the input so that the image quality can be compared relatively. Logarithmic transformation is applied to the raw pixel intensities of the detected image before reconstruction in the BP and SART methods. We assume a monoenergetic x-ray source and ignore the effects of scattering and beam hardening in this study, similar to the approach by Wu et al.2,24
We used four different phantoms to evaluate the reconstruction methods. The first phantom contains simulated microcalcification clusters and masses at different depths overlapping each other in a homogeneous background. This phantom is designed to evaluate image reconstruction quality for in-plane (X–Y) images and for off-plane (Z axis) blurring. As schematically shown in Fig. 2, the phantom consists of five 1-cm-thick breast shaped Lucite slabs. On top of the second, third, fourth, and fifth slabs, we placed different contrast objects including layers of thin circular aluminum foils to simulate masses (denoted as M*), calcium carbonate particles to simulate microcalcification clusters (denoted as C*), and high-contrast steel wires (denoted as W*). The objects on the same slab are distributed on a plane. The thicknesses of these objects are less than 0.3 mm, which is much thinner than the reconstruction slice thickness of 1 mm used in this study, and the objects can thus be considered to be located at the same depth. Two simulated masses (M1 and M3) and two groups of microcalcification clusters (C2 and C3), respectively, are overlapped along the Z-axis direction and separated by 1 cm.
The second phantom is the ACR mammography accreditation phantom which contains different simulated mammographic objects: nylon fibers, simulated calcification clusters, and simulated tumor-like masses. All objects are embedded in a wax block at about the same depth.
The third phantom was designed by us and custom-built by Computerized Imaging Reference Systems (CIRS), Inc., using breast-tissue-equivalent materials described in their website. It is composed of four 1-cm-thick breast-shaped slabs. The first and the fourth slabs are heterogeneous mixture of fibroglandular-tissue-mimicking material in a 80% fatty/20% glandular background. The second and the third slabs contain a homogeneous mixture of 50% fibroglandular and 50% fatty material. Several simulated spiculated masses are embedded in the second slab. On the upper surface of the third slab is a 6 × 5 array of drilled contrast-detail disk-shaped holes having diameters of 5, 4, 3, 2, 1, and 0.48 mm, and depths of 1, 0.8, 0.6, 0.4, and 0.2 mm. We will refer to this phantom as the breast C–D phantom in the following discussion.
A wire phantom is used to measure the relative sharpness of the reconstruction algorithms. A 220-µm-diameter steel wire in a thin plastic case is placed at about 2.5 cm above the surface of the breast support plate. Since the in-plane blurring artifacts of tomosynthesis occurs mainly along the tube-motion direction, the wire is oriented perpendicular to this direction. We will compare the normalized line profiles of the reconstructed wire images obtained with different methods.
To quantitatively evaluate the reconstructed image quality, we calculated the contrast, the root-mean-square (RMS) noise, and derived the CNR of selected features such as a mass or microcalcification at its in-focus plane. The CNR value is defined by
where Īfeature and ĪBG are the average pixel intensity of the feature and image background, respectively; and σBG is the RMS value of pixel intensity in the image background. The image background region for noise estimation is chosen as a 35 × 35-pixel region being far from all features and the boundaries of the imaged volume, and at the same slice as the feature under consideration. The average pixel intensity of a mass is calculated in a 35 × 35-pixel area enclosed within the relatively uniform central region of the mass, while that of a microcalcification in a 3 × 3-pixel area located approximately at the center of the microcalcification.
To evaluate image blur in the Z direction (perpendicular to the X–Y detector plane) of the reconstructed images and the artifact effect of features in the adjacent off-focus planes, we calculated an ASF, which is defined as the ratio of the CNR values between the off-focus layer and the feature layer:
where z0 is the slice location of the in-focus plane of the feature and z is the location of a plane of interest. The regions of interest (ROI) for analysis of mass and microcalcification, and the image background are the same as those described above for the calculation of CNR. A similar ASF has been used to describe the artifacts along the Z direction in a breast tomosynthesis study by the MGH research group,24 where the ASF was calculated as the ratio of contrasts only. Since the noise levels are expected to be different on the reconstructed images at different depths and the detectability of a feature is dependent more on the CNR than on the contrast alone, we choose to use the CNR rather than contrast alone in the ASF definition. The increase in noise with increasing depth of the image slice is caused by the cone-beam geometry of breast tomosynthesis which results in a decreasing x-ray intensity and fewer number of rays passing through a voxel as the depth increases. The reduced x-ray intensity increases quantum noise and the fewer number of rays makes the voxel updating more sensitive to the measurement noise.
To quantify the noise behavior of the reconstructed images as a function of spatial frequency, we calculate the relative NPS. NPS is an important measure to quantify the spatial-frequency-dependent noise level of x-ray imaging systems which are contributed by various noise sources such as the quantum noise, electronic noise, and Swank noise.43–46 In this study, we calculate one-dimensional (1D) NPS along both the X and Y directions on the reconstructed slices. The ROI for NPS determination is chosen as a 512 × 512-pixel homogeneous area on one tomosynthesized layer. Along each direction, we divide the region into multiple 512 × 16-pixel strips and the adjacent strips are overlapped by eight pixels, resulting in 63 samples. For each sample, the average pixel value is calculated along the sample to obtain a pixel intensity profile. No window function (implying a rectangular window) is used. To correct for the low frequency nonuniformities in the profile, a least-squares fit of a second-order polynomial is used to estimate the spatial profile of the pixel intensity and then the smoothed fit is subtracted from the original profile. A one-dimensional fast Fourier transform (FFT) is applied to the background-corrected noise profile. These fitting and FFT processes are performed separately for each sample. Finally, the 1D NPS is estimated by averaging the squared magnitude of the 1D FFT of all samples. The noise will depend on the exposure conditions and the phantom image chosen. However, since we are only interested in the relative noise performance from different reconstruction techniques in this study, we chose one phantom slice for this comparison and no normalization of the NPS was performed.
We define a line object spread function (LOSF) to compare the relative sharpness of high contrast linear objects in the tomosynthesized images using different reconstruction algorithms. The LOSF is conceptually similar to the line spread function (LSF) that is used to quantify the spatial resolution of an imaging system. However, since we do not use a very thin wire to simulate a delta function input and the wire image contains geometric unsharpness, the LOSF is different from an LSF. In this study, we are interested in measuring the relative blur of the reconstruction algorithms so the LOSF will be sufficient to serve the purpose.
For both the SART and ML-convex methods, we used the BP reconstruction results as initial values. In addition, a nonnegative constraint was applied to the reconstructions during the iterative process. For SART, three iterations were performed with decreasing step sizes of 0.3, 0.2, and 0.1 while a constant step size of 1.0 and up to eight iterations were evaluated for the ML-convex method. We determined these values experimentally to be representative ranges of step sizes and iterations that would provide reasonable reconstructed image quality in our applications.
ROIs from the reconstruction images containing the six simulated masses, four groups of microcalcification, and three steel wires, by the BP, SART, and ML-convex methods are shown in Fig. 3–Fig. 5, respectively. For all features, the SART images shown were obtained after the first iteration with a step size of 0.3 and the ML-convex images were obtained after seven iterations with a step size of 1.0. For comparison, the single-projection image of phantom 1 and the ROIs containing the overlapping masses M1 and M3, and microcalcification groups C2 and C3 are shown in Fig. 6.
The results showed that all three methods can reconstruct the features at their correct layers. The superimposed features along the Z direction, i.e., masses M1 and M3, and microcalcification clusters C2 and C3, were well separated. There were no obvious artifacts caused by overlapping structures, as can be seen by comparing the tomosynthesized slices of M1 and M3 in Fig. 3 with the single-projection image in Fig. 6(b), and those of C2 and C3 in Fig. 4 and Fig. 6(c). The effect of artifacts depends on the distance between the feature layers in the Z direction and the contrast of the features. In these examples, the feature layers are separated by 1 cm. Further analysis of the artifacts in the Z direction based on the ASF is presented below.
The BP method resulted in low-contrast features with smoothed edges while both the SART and ML-convex methods improved the conspicuity by enhancing the edge and the contrast. On the other hand, BP results produced smoothed image background with very low noise level, but both the SART and ML-convex methods, iterating from the BP results, strongly amplified the noise. This is due to the ill-posed nature of the inverse problem in tomosynthesis and no regularization has been explicitly applied to the algorithms. For high-attenuation features such as microcalcifications and wires, both iterative methods resulted in shadow regions extending from the object in the direction of the x-ray source motion.
To examine the effect of the iterative process on reconstruction image quality, we calculated the corresponding contrast, RMS noise and CNR values as the number of iterations increased. The phantom was imaged three times and the repeated measurements were averaged and their standard deviations estimated. The results for one of the masses, M5, and one microcalcification in group C3 (denoted by the same group name) are shown in Fig. 7. The values for the first three iterations of SART and the first eight iterations for ML-convex were included.
Figure 7 shows that the features reconstructed by the BP method are very low contrast but have relatively high CNR values because of the very low noise level. The SART method, with only one iteration, can significantly increase the contrast values, but simultaneously amplified the background noise to a high level, resulting in a low CNR value. Similar observations can be made for the ML-convex method, but the latter improved the contrast with a slower speed. The results in Fig. 7 suggested that SART with one iteration (λ = 0.3) can reach a similar CNR level to that of the ML-convex method with six or seven iterations. Similar conclusions can be obtained for both mass and microcalcification features. For this reason, we used the results from the SART method with one iteration (λ = 0.3) and the ML-convex method with seven iterations (λ = 1.0) in the following discussions.
Figure 8 shows the line profiles of both M5 and C3 for the BP, SART, and ML-convex methods. To get the line profile of the mass, five consecutive columns (161 pixels per column) were averaged with the central column passing through approximately the center of the mass. Similarly, for the microcalcification, three consecutive columns (21 pixels per column) were averaged with the central column passing through approximately the center of the microcalcification. The line profiles shown in Fig. 8 were mean removed and averaged over three repeated measurements. For both M5 and C3, the BP method resulted in a relatively smoother line profile while the SART and ML-convex methods gave sharper edges. This agrees with our observations from the reconstructed images. In addition, the SART and ML-convex methods gave very similar results for the selected features, as evident from their almost identical line profiles. Similar results were observed from analysis of other masses and microcalcifications although only one example of each is shown here.
Figure 9 shows the ASF curves of the selected mass and microcalcification. The layers with positive distance denote the image slices below the feature layer and vice versa. It is seen that the BP results have strong interplane blurring effect for the mass object, represented by a slowly decreasing ASF curve. Both the SART and ML-convex methods were superior in suppressing interplane blurring and the corresponding ASF curves dropped very quickly as the distance from the feature increased. For microcalcifications, all three methods have comparable ASF behaviors while BP gave a slightly better mean ASF curve.
The regular single-projection image of the ACR phantom is shown in Fig. 10(a). The tomosynthesized images of the ACR phantom feature layer (6 mm below the phantom surface) by the BP, SART, and ML-convex methods are shown in Figs. 10(b), 10(c), and 10(d), respectively. Similar to the results of phantom 1, the BP method produced images with low contrast features and low noise background. In the BP image, five fibrils, three groups of specks plus three specks in the fourth group, and three masses and part of the fourth mass may be seen on a good quality display monitor but it is difficult to see on the printed image. In the SART and ML-convex images, the same features are much more conspicuous with sharper edges and higher contrast. In addition, one more fibril and the whole fourth mass can be observed on both images while two more specks and the whole fourth group of specks can be observed for SART and ML-convex images, respectively.
We used the reconstructed ACR phantom images to calculate the 1D NPS along both X and Y directions. The tomosynthesized layer containing the homogeneous ROI for NPS estimation is 15 mm below the feature layer as shown in Fig. 10 (or 21 mm below the phantom surface). On this layer, the interplane artifacts from the feature layer were not observable. The choice of this phantom slice for the NPS comparison was somewhat arbitrary except that the slice should be far enough from the feature layer to avoid interplane artifacts. The measurement of the 1D NPS was repeated on three DTM samples of the ACR phantom at the same layer and the averages of the repeated measurement were compared.
The average 1D NPS for both X and Y directions in the same selected homogeneous area are shown in Fig. 11 for the three reconstruction methods. Since the digital detector has a pixel pitch of 0.1 mm, the Nyquist frequency occurs at 5 cycles/mm. The BP method produced much lower NPS level in the tomosynthesized slice, as also evident from the RMS noise estimate shown in Fig. 7. Both the SART (with one iteration) and the ML-convex (with seven iterations) methods significantly amplified the noise at all frequencies. The latter two have very similar NPS behaviors and the two curves in either direction are essentially indistinguishable.
The single-projection image of the breast C–D phantom is shown in Fig. 12(a). For comparison, we selected four reconstructed slices which are 6, 16, 20, and 34 mm below the breast C–D phantom surface that contain heterogeneous fibroglandular-tissue-mimicking material with different patterns (slices 1 and 4), several simulated spiculated masses on a homogeneous mixture (slice 2), and six columns of contrast-detail disk-shaped holes drilled in a homogeneous mixture (slice 3), respectively. The reconstruction results of the BP, SART, and ML-convex methods of the four slices are shown in Figs. 12(b)–12(d), 12(e)–12(g), 12(h)–12(j), and 12(k)–12(m), respectively. In the single-projection image, the contrast-detail disk array was obscured by the background structures and most of the disks were not visible except for a few large and high contrast ones. The simulated spiculated masses were almost totally invisible due to their low contrast and overlap with tissue structures. In the BP images, it is still very difficult to observe the speculated masses, while three rows of the drilled disks become visible. All disks on the first three rows except the smallest ones can be seen with strong windowing on the monitor. The homogeneous mixture background of these two feature slices (slices 2 and 3) were superimposed with severe streak artifacts along the x-ray source motion direction resulting from the tissue structure patterns both above and below. In contrast, both the SART and ML-convex methods can separate the overlapping tissue structures and the objects very well. It is much easier to observe the five spiculated masses and all disks in the first four rows including even the smallest ones in the first two rows. The objects have much sharper edges and higher contrast than those in the BP images. In slice 2, interplane artifacts from the higher contrast disks in the first row of the C–D array (slice 3) were visible. The background in both slices 2 and 3 also included some faint shadows from the heterogeneous tissue structures but the structured noise was much lower than the noise in the single-projection image and the BP images.
To compare the relative sharpness of high contrast line objects (edges) in the tomosynthesized images, we used the three methods to reconstruct images of a 220-µm-diameter steel wire that was placed about 2.5 cm above the surface of the breast support plate. The reconstructed layer in which the wire was in-focus was chosen as the feature layer. The wire was oriented perpendicular to x-ray source motion direction and the line profiles were taken in the direction perpendicular to the wire. Each line profile was derived from the average of four adjacent line profiles chosen from the same in-plane position for all three reconstruction methods. To determine the LOSF, we removed the background pixel intensity from the average line profile, and normalized its maximum intensity to 1. The wire phantom was imaged and the measurement of the LOSFs was repeated three times.
The resulting LOSF curves for one measurement are shown in Fig. 13. It is seen that the LOSFs of the SART and ML-convex methods had a similar level of broadening and both were narrower than that from the BP method. Similar trends of the LOSFs were observed for the repeated measurements. The mean and standard deviation of the full-width-athalf- maximum (FWHM) of the LOSFs for the BP, SART, and ML-convex methods estimated from the three measurements were 362±36, 342±27, and 279±19 µm, respectively. The undershoot lobes of the LOSFs for the SART and ML-convex methods are caused by the estimation of the average differences in the linear attenuation coefficients in the backprojection process which creates an effect similar to un-sharp masked filtering around edges. This may be one of the major factors contributing to the edge enhancement effects in the reconstructed DTM slices by the SART and ML-convex methods.
We have investigated three reconstruction algorithms, BP, SART, and ML-convex for DTMs using several phantoms. Preliminary results show that all methods can reconstruct the features in their correct layers and separate superimposed phantom structures along the Z direction. The two iterative methods chosen in this study, SART and ML-convex, are more effective in improving the conspicuity of object details and suppressing interplane blurring, as demonstrated in the reconstructed DTM images of the breast C–D phantom and by the improved ASF for masses measured with phantom 1. For the phantoms with homogeneous background as in phantoms 1 and 2, the BP method resulted in less noisy reconstruction and higher CNR values for low contrast objects than either the SART or the ML-convex method, but the latter methods provided stronger enhancement of the contrast of microcalcifications and edge sharpness of both masses and fibrils. Image noise for the SART and ML-convex methods increased as the number of iterations increased. With the breast C–D phantom, it was found that the BP method is inferior to the SART and ML-convex methods in terms of structured noise suppression and interplane blurring. Although the measured CNRs of low contrast objects for the BP method is higher than the CNRs for the two iterative methods using a phantom with a homogeneous background as shown in Fig. 7, the performance of the BP method in clinical images will be poor in the presence of structured background, as demonstrated in the breast C–D phantom images.
The selection of the relaxation parameter λ was experimentally determined in this study. To our knowledge, there is no general rule for choosing this parameter, and the most appropriate way is to determine it experimentally for a given application. For the ML-convex method, the use of constant step size, i.e., λ = 1.0, has been shown to provide satisfactory results in previous breast tomosynthesis reconstructions.2,47 In our study, we used the same strategy. For the SART method, the relaxation parameter (or step size) usually plays an important role. For all three phantom data investigated in this study, λ = 0.3 gave acceptable results after only one iteration. As the iterative process continued, the tomosynthesized slices quickly became too noisy. Thus, very few iterations in combination with adaptive relaxation parameters having fast decreasing values may be desirable in practical applications. The optimal strategy to choose relaxation parameters in SART for clinical cases remains to be investigated.
The computational burden of one iteration (all projection-view images have been processed exactly once) is the same order of magnitude for both SART and ML-convex methods, and most of the computation time is spent in the calculation of the simulated data (i.e., forward projection). The reason that the SART method has faster convergence speed than the ML-convex method can be attributed to, without considering the relaxation parameter, the update strategy of the unknowns. In the SART method, the linear attenuation coefficients are updated at each projection view; while for the ML-convex method, the update is performed after all views have been processed. The update information is over-smoothed in the ML-convex approach. It should be noted that some researchers have suggested using a block-iterative strategy in the ML-type methods to accelerate the convergence speed, which makes it close to the update strategy of SART.48,49 In this study, we chose the ordinary ML-convex implementation because it was used in previous investigations of breast tomosynthesis2,24 that might provide a reference for comparison, although a different imaging geometry and a different number of projection views were used in the previous studies.
The inverse problem of breast tomosynthesis reconstruction is an ill-posed problem due to the limited number of available projection views, the limited angular range of projection, and the large number of unknowns, which lead to serious amplification of noise in the data during the reconstruction. Regularization techniques may be useful for addressing this problem. We did not implement specific regularization methods. However, in our experiments, a degree of regularization was achieved through “early stopping,” i.e., the iteration was stopped before the solution became too noisy for both the SART and ML-convex methods, and through a non-negative constraint imposed on the unknown values. Further investigation on regularization methods may be needed to improve the image quality for breast tomosyn-thesis mammography.
The implementations of both the SART and ML-convex methods can be adapted to parallel computation algorithms. In both the forward projection and backprojection computation, the detector area can be divided into any reasonable number of subregions and the simulated measurements within each subregion can be calculated independently. The reconstruction of the entire volume can then be retrieved by stitching all the subregions together. This method has been used in breast tomosynthesis with the parallel ML-convex method,47 and chest tomosynthesis with algebraic reconstruction techniques.36
In this study, we have applied three representative reconstruction algorithms for the tomosynthesis of breast phantom data. Preliminary results show that for phantoms with homogeneous background, all methods can reconstruct the features in their correct layers and separate overlapping features. The BP method provided very smooth reconstructed images with low background noise, while the SART and ML-convex methods considerably enhanced the contrast and edges of the features but simultaneously amplified the image noise. The CNR values of simulated low-contrast objects for the BP method were higher than those for the SART and ML-convex methods. However, the two iterative methods can reduce the interplane blurring and artifacts with better ASF behaviors. For a contrast-detail phantom with heterogeneous tissue-mimicking background, the BP method had blurring artifacts in the x-ray source motion direction that obscured the contrast-detail objects, while the other two methods can significantly improve object conspicuity by removing the overlapping structures. With a properly selected relaxation parameter, the SART method with one iteration can provide tomosynthesized images comparable to those obtained by using the ML-convex method with seven iterations, when the BP results were used as initialization in both methods. Future work will be conducted to evaluate the effects of various parameters including initialization, step size, number of projection views, angular increments and ranges, and regularization on the quality of DTM images. In addition, studies will be performed to compare the image quality of patient DTM images using different reconstruction methods.
This work is supported by USPHS Grant Nos. CA91713 and CA95153, and U. S. Army Medical Research and Materiel Command Grant No. DAMD 17-02-1-0214. The digital tomosynthesis system was developed by the GE Global Research Group through the Biomedical Research Partnership (USPHS Grant No. CA91713) collaboration.