Use of the Hotelling observer to optimize image reconstruction in digital breast tomosynthesis

Abstract. We propose an implementation of the Hotelling observer that can be applied to the optimization of linear image reconstruction algorithms in digital breast tomosynthesis. The method is based on considering information within a specific region of interest, and it is applied to the optimization of algorithms for detectability of microcalcifications. Several linear algorithms are considered: simple back-projection, filtered back-projection, back-projection filtration, and Λ-tomography. The optimized algorithms are then evaluated through the reconstruction of phantom data. The method appears robust across algorithms and parameters and leads to the generation of algorithm implementations which subjectively appear optimized for the task of interest.


Introduction
In recent years, digital breast tomosynthesis (DBT) has continued to gain attention as a promising modality for breast imaging. By combining projection images acquired at varying angles, the effects of tissue superposition can be ameliorated, leading to improved visualization of masses and a decrease in false positives for mass detection. [1][2][3][4] Since the angular range used in tomosynthesis is small, DBT acquisitions contain insufficient data to allow conventional tomographic image reconstruction methods to yield fully three-dimensional (3-D) images. Instead, a range of variants of analytic reconstruction methods from x-ray CT, such as filtered back-projection (FBP), has been developed for DBT in order to obtain quasi-3-D images. On the one hand, the inability to achieve full object representation in 3-D means that without incorporating prior object information, it is not possible to obtain an image which uniquely corresponds to the projection data. On the other hand, this nonuniqueness implies a certain freedom in designing an image reconstruction algorithm, with even simple back-projection yielding images of potential utility. While there have been some heuristics to guide DBT reconstruction research, to date DBT algorithm design has been largely guided by assessment of image artifacts. 1,4 In this work, we present a method for objectively optimizing analytic reconstruction algorithms in DBT. The method is based on the Hotelling observer (HO) 5 and optimizes a given algorithm for performance of a signal detection task. We focus on two tasks in this work: the detection of microcalcifications, small deposits of calcium which can indicate malignancy, and the detection of small low-contrast disks, intended to mimic subtle lesions. Previous work by others has investigated the use of the HO for optimization of DBT acquisition parameters, [6][7][8][9] benchmarking DBT acquisition relative to mammography, 8,10 and, similar to this work, exploring the impact of image reconstruction. 6,11 Often, the HO's performance for a given task is estimated using a collection of sample images, resulting in statistically variable estimates of the HO's figure of merit. Since we wish to perform optimization of several reconstruction parameters, for efficiency, we construct an approximation of HO performance which does not rely on samples of noisy images, and is therefore nonstochastic. Most previous studies in constructing nonstochastic HO performance estimates rely on channels, 12 such as Laguerre-Gauss channels, 13 which rely on assumptions of symmetry in the signal of interest in order to reduce the dimensionality of the HO metric computation. Instead, in this work the HO signal-to-noise ratio (SNR) is obtained in the spatial domain without the use of channels. Rather, extending previous work, [14][15][16] the dimensionality of the HO system of equations is reduced by considering microcalcification detection only within a restricted region of interest (ROI). Meanwhile, the effect of reconstruction on the signal and quantum noise is captured by explicitly modeling the effect of each linear operation on a Gaussian detector noise model. We demonstrate that the proposed method can enable efficient simultaneous optimization of multiple reconstruction parameters. This is potentially useful for DBT system development, where the interplay between various reconstruction parameters and their collective effect on image quality can be difficult to predict a priori. This method, in combination with more conventional assessment methods (e.g., artifact evaluation), can provide a fuller picture of the utility of a given algorithm implementation for DBT.
Here, we demonstrate the use of the HO to optimize four reconstruction algorithms, three of which are variants of the conventional FBP algorithm, with the fourth algorithm being simple back-projection without filtration. These algorithms each contain several parameters which can be difficult to determine manually. After optimizing the algorithms for two signal detection tasks, we subjectively evaluate the implementations we obtain using phantom data from a Hologic tomosynthesis unit. The quality of these reconstructed images is assessed in light of the values obtained for the HO SNR figure of merit for each reconstruction algorithm.
Since our primary motivation in this work is to facilitate the development and implementation of reconstruction algorithms, we place a particular emphasis on designing an assessment method which is simple and efficient, so that extensive parameter sweeps can be performed to yield optimized reconstruction algorithms. Once the large parameter space which exists for reconstruction algorithms is reduced by this simplified model, more realistic modeling could be used to evaluate the entire imaging system and compare relatively few algorithms or parameter choices.

Digital Breast Tomosynthesis System Model
We will use the convention that two-dimensional (2-D) or 3-D images are represented by one-dimensional (1-D) vectors, indexed by a single variable, which can be obtained by lexicographical ordering. We begin the development of our DBT system model by operating in the projection data domain after the negative logarithm has been applied. The i'th element of the projection data g can then be modeled as a line integral through the continuous object f E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 4 7 4 g i ¼ where the domain of integration L is defined as the intersection of the i'th ray with the compact support of the object, s i denotes the source position, andθ i is a unit vector in the direction of the detector element corresponding to the i'th ray. Further, we assume that the noise in each detector reading can be modeled as independent, additive Gaussian noise via the stochastic term n i . In order to model the effect of finite detector pixels, the model in Eq. (1) is applied with 16-fold subsampling within each detector element. Finally, a 0.4-mm focal spot size is modeled by convolving the projection data with a 0.4-mm rect function, scaled to account for geometric magnification. For simplicity, we will adopt a noise model which is strictly Gaussian and assumes uncorrelated noise in the detector elements. Under this assumption, the data covariance matrix, denoted K g and given by K g ¼ E½ðg −ḡÞðg −ḡÞ T is diagonal, consisting only of variances with all covariances between detector elements equal to 0. Here and elsewhere, variables with bars denote averages and the superscript T denotes the transpose operator. The specific noise model we consider is taken from Ref. 17. In particular, the noise model begins with the assumption that the variance in the postlogarithm data goes as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 2 0 0 Varfg i g ¼ whereN 0 is a constant denoting the average number of photons at the i'th ray in a blank scan, andN i is the mean number of photons from the i'th ray that are transmitted through the object and reach the detector. Relating the termN i to the mean line integralḡ i gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 1 0 1 ðK g Þ i;j ¼ We have setN 0 ¼ 10 5 in this work; however, this choice does not impact the optimal reconstruction parameter settings, as a change inN 0 will scale all of the calculated HO SNR values by the same multiplicative factor. Electronic noise, detector cross talk, and detector efficiency were not modeled in this work.
In order to simplify the observer model and facilitate efficient computation while still capturing the essential imaging model, we only consider reconstruction within a specific ROI. This ROI is defined as the intersection of two planes. The first plane is parallel to the detector and contains the signal of interest. We restrict the observer to this plane because we make the assumption that only a single 2-D image is used for performing classification. This image can be either a single slice or as in the case of the back-projection filtration (BPF) algorithm a weighted average of several slices. The second plane is perpendicular to the detector, contains the signal, and is parallel to the scan trajectory. This plane is selected because it is likely to contain the majority of the signal information, since the signal is back-projected within this plane. The final ROI is then a 1-D profile through the signal of interest, running the full width of the reconstructed image in the scan direction. Once parameters such as pixel size are obtained along the scan trajectory direction (y-direction) using the proposed HO method, these are then matched within each slice in the direction from the chest wall to the nipple (x-direction).
In order to evaluate the performance of the optimized reconstruction algorithms, projection data of a Gammex 156 mammography phantom acquired using a Hologic tomosynthesis unit were used. Therefore, the configuration for the simulation portion of the study mimicked that system, with a source-toimage distance of 70 cm and 0.14 mm detector pixels. There were 15 total simulated projections, spaced in 1 deg increments, and for the simulations the axis of rotation was set to be even with the detector plane. Microcalcifications of 0.16 mm were simulated as Gaussians with maximum values equal to the mean linear attenuation coefficient of calcium for an 80-kVp tube spectrum. The smallest low-contrast mass in the Gammex phantom was modeled as a 0.25-cm uniform disk with 5% contrast relative to the background. The thickness of the disk was set to 0.25 cm, resulting in a contrast of 0.57% relative to background in the projection data. The simulated background was uniform, corresponding to the background of the Gammex phantom. As described in Sec. 1, the purpose of investigating reconstructions of real data is to perform subjective visual validation of the simplified HO model in light of the physical factors which are present in real data, but are not included in the simulation methodology.

Image Reconstruction Algorithms
In addition to pixel-driven back-projection, three reconstruction algorithms, all based on FBP, are considered in this work. The first algorithm is most similar to an FBP implementation commonly used in x-ray CT, with application of a ramp filter and Hanning apodization window applied before back-projection. The second algorithm we optimize is related to FBP, but performs the filtering step after back-projection. We will refer to this implementation as BPF to distinguish it from FBP; however, elsewhere this approach is also referred to as FBP. 18 Finally, we will also implement and optimize a reconstruction method known as Λ-tomography 19 or local tomography, 20,21 which has been proposed as a means of performing local image reconstruction in tomosynthesis.
Each of these four reconstruction algorithms is linear and discrete-to-discrete since it transforms the discrete projection data to discrete image vectors. The algorithm can then be thought of as a matrix, A, which produces an image vector y via the transform y ¼ Ag. Further, it will be beneficial for our purposes to decompose each algorithm into each of its constituent linear processing steps, so that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; 6 5 where there are N total processing steps in the algorithm, such as discrete Fourier transforms (DFTs), convolutions, etc. We can then obtain the useful result that the image covariance matrix K y is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; where the superscript † denotes the conjugate transpose and the limits of the first product imply reverse ordering of the product elements, i.e., A † N A † N−1 : : : A † 1 . Given implementations of each algorithm step A i , it is typically straightforward to construct an implementation of A † i . While K y is typically too large to store directly in computer memory, Eq. (5), together with the data noise model of Eq. (3), provides a convenient linear operator implementation for computing inner products with the image covariance matrix. This will provide the basis for iterative computation of the HO figure of merit described in Sec. 2.3.

Filtered back-projection algorithm
Here, we briefly summarize the fan-beam FBP algorithm implementation used in this work. Our description is based on that contained in Ref. 22. Fan-beam FBP reconstruction follows a three-step process involving weighting, filtration, and backprojection. First, each g i is multiplied by a weighting factor. This can be written in matrix form as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 2 9 8 g 0 ¼ Wg; where W is a diagonal weighting matrix whose elements are given by wðnÞ ¼ D∕ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Here, D is the source-todetector distance, n is an integer denoting the detector element where n ¼ 0 corresponds to the central ray, and d u is the detector pixel size. The matrix W is diagonal because each of the weights for the various detector pixels is applied independently.
Next, the modified projection data vector g 0 is filtered with a ramp filter, which can be implemented in the DFT domain by multiplication with another diagonal matrix, R, containing the DFT-domain ramp kernel. Likewise, apodization filters (such as a Hanning filter in this work) can be applied as another diagonal matrix in the DFT-domain. Finally, one performs a weighted back-projection of each of the filtered projections onto a discrete image array via a back-projection matrix operator B, resulting in the reconstructed image pixels which are elements of y. The full reconstruction operation can then be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 8 9y where H is a diagonal matrix with the Hanning window in the diagonal elements, R and W are again the ramp-filter and weighting matrices, and B and F 1-D are matrices representing back-projection and the discrete FFT operation, respectively. The subscript 1-D highlights the fact that the filtering is performed in the projection data domain only along the scan direction. The Hanning window is in turn given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 3 2 6 ; 6 7 5 where ν u is the detector spatial frequency, ν c is a parameter which controls the strength of filtering, and d u is again the detector pixel size. We will refer to ν c as the spectral filter cutoff, which is normalized so that ν c ¼ 1 corresponds to a cutoff at the Nyquist frequency of the detector. As pointed out by Ref. 18, the sampling properties of tomosynthesis introduce artifacts that are far more significant than those which result from neglecting the exact fan-beam weighting. Therefore, for simplicity, we will make a final modification of the FBP algorithm and set W equal to the identity matrix and perform unweighted pixeldriven back-projection in our simulations. For FBP, then, there are three significant parameters which require optimization: the in-plane pixel size, the slice thickness, and ν c , the cutoff for the Hanning filter.

Back-projection filtration algorithm
As an alternative to the more CT-like FBP algorithm just presented, Mertelmeier 18 describes a method for performing all of the filtering necessary for reconstruction in the image domain after back-projection. We refer to this as the BPF algorithm. After applying the back-projection matrix B to the projection data vector g, a 2-D DFT, denoted here by the matrix F 2-D , is applied in the image domain. The image spatial frequency version of the ramp filter is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 3 5 9 R BPF ≈ 2αjν y j; where ν y is the spatial frequency along the source trajectory direction (y-direction), and α is known as the tomo-angle in a scan with an angular range from −α to α. As in the FBP case, smoothing can be implemented via a Hanning filter; however, while the detector-domain smoothing affected the image along the scan direction (y-direction) and the slice direction (z-direction), image-domain filtering can be independently performed in the yand zdirections. (As with pixel-size, the smoothing applied to the x-dimension is set to match the y-dimension to preserve image appearance within a given slice.) For the BPF algorithm, we will refer to the scale factor in Eq. (8) as the spectral filter width or the slice filter width for the yand z-dimensions, respectively. The final reconstructed image is then given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 1 7 3 where H spectrum and H slice are given by Eq. (8) with ν y and ν z , respectively, replacing ν u . The matrices B and F 2-D again denote back-projection and the Fourier transform, and R BPF is defined in Eq. (9). For the BPF algorithm, as with FBP, pixel size and slice thickness are free parameters that must be determined, as are the spectral and slice filter cutoff frequencies.

Λ-tomography
The final algorithm we will investigate is Λ-tomography, [19][20][21] which reconstructs an image corresponding to a Laplacianfiltered version of the object, rather than a representation of the object itself. Λ-tomography has been advocated for tomosynthesis due to its ability to lessen artifacts due to projection truncation while preserving edges and regions of object uniformity. 19 While several variants exist, in this study we implement the most basic form, which simply replaces the DFT-domain ramp filter with a negative numerical second derivative, where the second derivative of the i'th detector element u i is approximated as u iþ1 − 2u i þ u i−1 . For the first and last detector elements (i ¼ 1 and i ¼ i max ), we assume that u 0 ¼ u 1 and u i max þ1 ¼ u i max . Denoting the matrix form of this operator as L, we obtain the following for the Λ-tomography algorithm E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; where, similar to FBP, we have again allowed for spectral filtering in the data domain via the Hanning filter matrix H. The matrices L, F 1-D , and B are again the negative numerical second derivative, the FFT, and back-projection, respectively. The parameters to be optimized for Λ-tomography are essentially the same as those for FBP: slice thickness, pixel size, and Hanning filter cutoff.

Hotelling Observer Signal-to-Noise Ratio Metric
The HO's figure of merit is the HO SNR, 5 given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 4 3 1 SNR 2 ¼ w T Δȳ; (12) where for zero-mean additive noise, Δȳ is equivalent to a noisefree reconstruction of the signal of interest and w is the Hotelling template, defined implicitly through the equation E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 3 6 6 K y w ¼ Δȳ: (13) Recall that the term K y is the image covariance matrix, given by Eq. (5). Typically, some dimension-reducing operator is appended to the reconstruction algorithm so that the matrix K y can be efficiently computed and inverted. Here, in lieu of conventional channels, we reduce the dimensionality of K y by restricting the reconstruction to an ROI, as described in Sec. 2.1. The resulting image covariance is still too large to store directly in computer memory, so Eq. (13) is instead solved iteratively, using the method of conjugate gradients. For the present work, the solution via conjugate gradients is obtained with sufficient efficiency that a grid search of the algorithm parameter space is used to obtain optimal parameters for microcalcification detection. As a check, the HO SNR in the projection data domain is also computed, as this is an upper bound on the SNR estimated by our postreconstruction HO implementation. The squared ratio of the image domain HO SNR to the projection-domain HO SNR is then reported as an efficiency value.
Of key importance to the present work is that, as constructed, the SNR of Eq. (12) is not defined in terms of estimates derived from sample images. By contrast, the HO SNR we obtain is nonstochastic. In most studies of HO performance, HO SNR is estimated from samples, so that bias and variance of the estimates must be considered when comparing our approach with more conventional methods for computing HO SNR. For a more in-depth discussion of this issue and experimental comparisons of nonstochastic ROI-based methods and stochastic samplebased methods, see Ref. 14.

Results
The results for optimal efficiency values and parameter settings for microcalcification detection and low-contrast signal detection are shown in Tables 1 and 2, respectively, for each of the algorithms considered. Interestingly, despite restriction to an ROI, most algorithm implementations were extremely efficient with regard to preserving signal detectability from the projection data domain, with over 99% of HO SNR being preserved postreconstruction in the case of simple back-projection. The HO efficiency was found to vary with slice thickness primarily because of the resulting location of the slice center relative to the microcalcification signal; therefore, we also allowed the slice thickness to be a free parameter for optimization by the HO. For cases where there seemed to be no SNR penalty for decreasing pixel size, the largest pixels with near-optimal SNR were selected. Figure 1 illustrates trends of HO efficiency with varying pixel size and filter width for BPF, FBP, and Λ-tomography. For the low-contrast disk detection, HO efficiency was insensitive to pixel size across a wide range. Drastic jumps in efficiency of microcalcification detection with varying pixel size are an effect of discretization, with performance increasing and decreasing as the modeled signal moves in and out of pixel centers. Similar trends were seen for slice thickness (not shown). Note that in many cases, the HO efficiency is relatively insensitive with a wide plateau near an optimal parameter setting. It is also interesting to observe the trade-off in optimizing reconstruction parameters for each of the two tasks. For instance, performance trends for each task are opposite with respect to filter cutoff. For the low-contrast detection task, heavier filtering does not substantially degrade the contrast, but reduces noise magnitude, improving performance. Meanwhile, low-frequency filter cutoffs strongly degrade microcalcification detection performance, as these signals have substantial high spatial frequency content. Intersections in each of these plots define compromise parameter settings for algorithm implementations that are equally efficient at preserving information for each task. In addition to phantom images reconstructed using algorithms optimized for a single task, we next present images reconstructed with parameter settings of equal HO efficiency for each task. Figure 2 shows reconstructed images of the 0.32 mm group of simulated microcalcifications (specks) in the Gammex phantom. This is the smallest group for any of the reconstruction methods where all six specks are visible. The left side of the figure corresponds to the reconstruction implementations deemed optimal by the HO for microcalcification detection, while for the images on the right side, the spectral and slice filter widths were set for equal HO efficiency for the two tasks (nothing was changed for simple back-projection). The correspondence between images and algorithms is given in the figure caption. The interpolation artifacts in the FBP images are a result of pixel size being mismatched with the detector pixel size; however, the HO is insensitive to these artifacts. Likewise, Fig. 3 shows similar reconstruction results for the central signal in the group of 0.24 mm specks. Since the algorithms result in very different gray values, a dynamic windowing scheme is used to display the images.
The images in Fig. 4 show similar results, but for algorithms optimized for the low-contrast disk detection task. Again, since optimal pixel size and slice width were insensitive to the choice  of task, nothing is changed for simple back-projection between the two subfigures. The disk-optimized Λ-tomography implementation demonstrates significant ringing artifacts near object edges from over-filtering, but the HO is again insensitive to these artifacts. By inspection of Fig. 4, it is apparent that the low-contrast inserts in the Gammex phantom have a wide range of appearances depending on the reconstruction algorithm used.
In general, although the images obtained with different algorithms have very different appearances and parameters that vary substantially (e.g., pixel size between FBP and Λ-tomography), the reconstructions of the phantom data appear reasonable in light of the quantitative results obtained for HO efficiency in each case. Since, based on the efficiency values in the two tables, the HO was able to select algorithm implementations which were nearly equivalent in the context of a specific task, the method demonstrated here could be used as a means of more holistic algorithm assessment based on artifacts and other properties of images like those shown in Figs. 2-4, while keeping signal detectability fixed.

Conclusions
In this work, we have proposed an implementation of the HO which operates in the spatial domain without the use of channels in order to optimize linear reconstruction algorithms in DBT. Instead of using channels, only a single dimension of the image volume is utilized for HO assessment in order to reduce dimensionality, and the Hotelling template is subsequently obtained iteratively. Our results demonstrate that up to roughly 99% of the information in the projection data domain which is useful to the HO is still contained in the 1-D domain used by our HO implementation. Further, subjective evaluation of the algorithm implementations obtained with this method suggest that, while the HO is insensitive to some artifacts, this method leads to images of quality which is at least comparable to that obtained through time-consuming manual parameter tuning.
Some limitations of the present work are that, since the methodology proposed was designed to be as simple as possible for efficient algorithm optimization, many choices in modeling the system and task lack realism. This suggests that, as it stands, the approach demonstrated here is not adequate for end-point evaluation of an imaging system for the purpose of determining performance for tasks with immediate clinical relevance. Instead, the methods we put forward are more suited to component-level optimization of the reconstruction algorithm in the development stage of system design. More realism may be possible for the lightweight HO implementation we propose, such as nearest neighbor detector pixel correlation, different signal models, overlapping structures or tumors, or possibly a simple random background model. The main purpose of the lightweight HO is to reduce the large parameter space which exists for reconstruction algorithms. Realistic anthropomorphic model observers, which are more costly to evaluate, could then be used to select between comparatively few parameter settings or algorithms.