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

Adrian A. Sánchez; Emil Y. Sidky; Xiaochuan Pan

doi:10.1117/1.JMI.3.1.011008

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

Adrian A. Sánchez; Emil Y. Sidky; Xiaochuan Pan

Author Affiliations +

J. of Medical Imaging, 3(1), 011008 (2015). doi:10.1117/1.JMI.3.1.011008

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.

Sánchez, Sidky, and Pan: Use of the Hotelling observer to optimize image reconstruction in digital breast tomosynthesis

1. 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.¹2.3.^–⁴ 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.¹^,⁴

In this work, we present a method for objectively optimizing analytic reconstruction algorithms in DBT. The method is based on the Hotelling observer (HO)⁵ 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,⁶7.8.^–⁹ benchmarking DBT acquisition relative to mammography,⁸^,¹⁰ and, similar to this work, exploring the impact of image reconstruction.⁶^,¹¹

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,¹² such as Laguerre–Gauss channels,¹³ 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,¹⁴15.^–¹⁶ 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.

2. Methods

2.1.

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$

(1)

g_{i} = \int_{L} f (s_{i} + l {\hat{θ}}_{i}) d l + n_{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

{\hat{θ}}_{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 - \bar{g}) {(g - \bar{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

(2)

Var {g_{i}} = \frac{1}{{\bar{N}}_{i}} + \frac{1}{{\bar{N}}_{0}},

where

{\bar{N}}_{0}

is a constant denoting the average number of photons at the

i

’th ray in a blank scan, and

{\overline{N}}_{i}

is the mean number of photons from the

i

’th ray that are transmitted through the object and reach the detector. Relating the term

{\bar{N}}_{i}

to the mean line integral

{\bar{g}}_{i}

gives

(3)

{(K_{g})}_{i, j} = {\begin{matrix} \frac{e^{{\bar{g}}_{i}} + 1}{{\bar{N}}_{0}} & : i = j \\ 0 & : else . \end{matrix}

We have set

{\bar{N}}_{0} = 10^{5}

in this work; however, this choice does not impact the optimal reconstruction parameter settings, as a change in

{\bar{N}}_{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-to-image 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.

2.2.

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.¹⁸ Finally, we will also implement and optimize a reconstruction method known as $Λ$ -tomography¹⁹ or local tomography,²⁰^,²¹ 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 = A g$ . Further, it will be beneficial for our purposes to decompose each algorithm into each of its constituent linear processing steps, so that

(4)

A = \prod_{i = 1}^{N} A_{i},

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

(5)

K_{y} = (\prod_{i = N}^{1} A_{i}^{†}) K_{g} (\prod_{i = 1}^{N} A_{i}),

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}^{†} \dots 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.

2.2.1.

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 back-projection. First, each $g_{i}$ is multiplied by a weighting factor. This can be written in matrix form as

(6)

g^{'} = W g,

where

W

is a diagonal weighting matrix whose elements are given by

w (n) = D / \sqrt{D^{2} + n^{2} d_{u}^{2}}

. Here,

D

is the source-to-detector 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^{'}$ 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

(7)

y = A g = B F_{1 ‐ D}^{- 1} H R F_{1 ‐ D} W g,

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

(8)

H (ν_{u}) = 0.5 [1 + \cos (\frac{2 d_{u} ν_{u}}{ν_{c}})],

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 pixel-driven 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.

2.2.2.

Back-projection filtration algorithm

As an alternative to the more CT-like FBP algorithm just presented, Mertelmeier¹⁸ 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

(9)

R_{BPF} \approx 2 α | ν_{y} |,

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

y

- and

z

- directions. (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

y

- and

z

-dimensions, respectively. The final reconstructed image is then given by

(10)

y = F_{2 ‐ D}^{- 1} H_{spectrum} H_{slice} R_{BPF} F_{2 ‐ D} B g,

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.

2.2.3.

$Λ$ -tomography

The final algorithm we will investigate is $Λ$ -tomography,¹⁹20.^–²¹ which reconstructs an image corresponding to a Laplacian-filtered 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.¹⁹ 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} - 2 u_{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

(11)

y = A g = B F_{1 ‐ D}^{- 1} H F_{1 ‐ D} L g,

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.

2.3.

Hotelling Observer Signal-to-Noise Ratio Metric

The HO’s figure of merit is the HO SNR,⁵ given by

(12)

{SNR}^{2} = w^{T} Δ \bar{y},

where for zero-mean additive noise,

Δ \bar{y}

is equivalent to a noise-free reconstruction of the signal of interest and

w

is the Hotelling template, defined implicitly through the equation

(13)

K_{y} w = Δ \bar{y} .

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 sample-based methods, see Ref. 14.

3. 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.

Table 1

Optimal parameters for microcalcification detection for the four reconstruction methods considered in this work, as determined by the HO implementation we propose. Also presented are the efficiency values of the HO in the image domain relative to the HO SNR operating directly on the projection data.

	Efficiency (%)	Optimal slice thickness (mm)	Optimal y-pixel size (mm)	Spectral filter cutoff	Slice filter cutoff
Back-projection	99.70	1.125	0.085	N/A	N/A
BPF	96.89	1.184	0.0854	1.28	1.68
FBP	97.06	1.41	0.132	1.2	N/A
Λ-tomography	87.25	1.184	0.04457	1.13	N/A

Table 2

Same as Table 1, but for the low-contrast detection task.

	Efficiency (%)	Optimal slice thickness (mm)	Optimal y-pixel size (mm)	Spectral filter cutoff	Slice filter cutoff
Back-projection	99.85	1.125	0.085	N/A	N/A
BPF	91.03	1.184	0.178	0.108	0.573
FBP	92.31	1.41	0.178	0.121	N/A
Λ-tomography	81.19	1.184	0.421	0.0393	N/A

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.

Fig. 1

Parameter sweeps of pixel size, spectral filter, and slice filter (if applicable) for three algorithms investigated—BPF, FBP, and $Λ$ -tomography. Results are shown for both the microcalcification detection task and the detection of the 0.25 cm low-contrast phantom insert. As a given parameter varies, all other parameters are kept fixed at their optimum values for the respective 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.

Fig. 2

Clusters of 0.32 mm specks from the Gammex phantom, obtained using (a) algorithms with parameters optimized by the HO for microcalcification detection and (b) algorithms optimized for equal efficiency for both tasks. For each subfigure, the images correspond to (top left) back-projection, (top right) FBP, (bottom left) BPF, and (bottom right) $Λ$ -tomography.

Fig. 3

Same as Fig. 2, but only the central speck from the 0.24 mm group of specks is shown. The 0.24 mm speck is located in the center of each image. As before, for each subfigure, the images correspond to (top left) back-projection, (top right) FBP, (bottom left) BPF, and (bottom right) $Λ$ -tomography.

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.

Fig. 4

(a) Reconstructed images of the full Gammex phantom corresponding to algorithms optimized for detection of the 0.25 cm low-contrast disk and (b) algorithms optimized for equal HO efficiency between the two tasks. Again, the arrangement of the subfigures corresponds to (top left) back-projection, (top right) FBP, (bottom left) BPF, and (bottom right) $Λ$ -tomography.

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.

4. 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.

Acknowledgments

The authors would like to thank Hologic for providing the phantom data. This work was supported in part by NIH R01 Grant Nos. CA158446, CA182264, and EB018102. The contents of this article are solely the responsibility of the authors and do not necessarily represent the official views of the National Institutes of Health.

References

1.

III J. T. Dobbins and D. J. Godfrey, “Digital x-ray tomosynthesis: current state of the art and clinical potential,” Phys. Med. Biol. 48, R65 (2003).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/48/19/R01Google Scholar

2.

III J. T. Dobbins, “Tomosynthesis imaging: at a translational crossroads,” Med. Phys. 36, 1956–1967 (2009).MPHYA60094-2405http://dx.doi.org/10.1118/1.3120285Google Scholar

3.

I. Sechopoulos, “A review of breast tomosynthesis. Part I. The image acquisition process,” Med. Phys. 40(1), 014301 (2013).MPHYA60094-2405http://dx.doi.org/10.1118/1.4770279Google Scholar

4.

I. Sechopoulos, “A review of breast tomosynthesis. Part II. Image reconstruction, processing and analysis, and advanced applications,” Med. Phys. 40(1), 014302 (2013).MPHYA60094-2405http://dx.doi.org/10.1118/1.4770281Google Scholar

5.

H. H. Barrett and K. J. Myers, Foundations of Image Science, John Wiley & Sons, Inc., Hoboken, New Jersey (2004).Google Scholar

6.

I. Reiser, B. A. Lau and R. M. Nishikawa, “Effect of scan angle and reconstruction algorithm on model observer performance in tomosynthesis,” in Digital Mammography, pp. 606–611, Springer, Heidelberg (2008).Google Scholar

7.

A. S. Chawla et al., “A mathematical model platform for optimizing a multiprojection breast imaging system,” Med. Phys. 35(4), 1337–1345 (2008).MPHYA60094-2405http://dx.doi.org/10.1118/1.2885367Google Scholar

8.

I. Reiser and R. M. Nishikawa, “Task-based assessment of breast tomosynthesis: effect of acquisition parameters and quantum noise,” Med. Phys. 37(4), 1591–1600 (2010).MPHYA60094-2405http://dx.doi.org/10.1118/1.3357288Google Scholar

9.

S. Young et al., “A virtual trial framework for quantifying the detectability of masses in breast tomosynthesis projection data,” Med. Phys. 40(5), 051914 (2013).MPHYA60094-2405http://dx.doi.org/10.1118/1.4800501Google Scholar

10.

S. Park et al., “A statistical, task-based evaluation method for three-dimensional x-ray breast imaging systems using variable-background phantoms,” Med. Phys. 37(12), 6253–6270 (2010).MPHYA60094-2405http://dx.doi.org/10.1118/1.3488910Google Scholar

11.

R. Zeng et al., “Evaluating the sensitivity of the optimization of acquisition geometry to the choice of reconstruction algorithm in digital breast tomosynthesis through a simulation study,” Phys. Med. Biol. 60(3), 1259 (2015).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/60/3/1259Google Scholar

12.

K. J. Myers and H. H. Barrett, “Addition of a channel mechanism to the ideal-observer model,” J. Opt. Soc. Am. A 4(12), 2447–2457 (1987).http://dx.doi.org/10.1364/JOSAA.4.002447Google Scholar

13.

B. D. Gallas and H. H. Barrett, “Validating the use of channels to estimate the ideal linear observer,” J. Opt. Soc. Am. A 20(9), 1725–1738 (2003).http://dx.doi.org/10.1364/JOSAA.20.001725Google Scholar

14.

A. A. Sanchez, E. Y. Sidky and X. Pan, “Region of interest based Hotelling observer for computed tomography with comparison to alternative methods,” J. Med. Imaging 1(3), 031010 (2014).JMEIET0920-5497http://dx.doi.org/10.1117/1.JMI.1.3.031010Google Scholar

15.

A. A. Sanchez, E. Y. Sidky and X. Pan, “Task-based optimization of dedicated breast CT via Hotelling observer metrics,” Med. Phys. 41(10), 101917 (2014).MPHYA60094-2405http://dx.doi.org/10.1118/1.4896099Google Scholar

16.

A. A. Sanchez, “Task-based optimization of computed tomography imaging systems,” Doctoral Thesis, University of Chicago, Chicago, Illinois (2015).Google Scholar

17.

H. H. Barrett and W. Swindell, Radiological Imaging: The Theory of Image Formation, Detection, and Processing, Academic Press, San Diego (1996).Google Scholar

18.

T. Mertelmeier, “Filtered backprojection-based methods for tomosynthesis image reconstruction,” in Tomosynthesis Imaging, , I. Reiser, S. Glick and W. R. Hendee, Eds., Imaging in Medical Diagnosis and Therapy, pp. 101–106, Taylor & Francis Group, Boca Raton (2014).Google Scholar

19.

J. Frikel and E. T. Quinto, “Characterization and reduction of artifacts in limited angle tomography,” Inverse Prob. 29(12), 125007 (2013).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/29/12/125007Google Scholar

20.

A. Faridani, E. L. Ritman and K. T. Smith, “Local tomography,” SIAM J. Appl. Math. 52(2), 459–484 (1992).SMJMAP0036-1399http://dx.doi.org/10.1137/0152026Google Scholar

21.

A. Faridani et al., “Local tomography II,” SIAM J. Appl. Math. 57(4), 1095–1127 (1997).SMJMAP0036-1399http://dx.doi.org/10.1137/S0036139995286357Google Scholar

22.

A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York (1988).Google Scholar

Biography

Adrian A. Sánchez received his BS degree in physics from the University of Tennessee, Knoxville, in 2010, and his PhD in medical physics from the University of Chicago in 2015. He is currently a diagnostic imaging physics resident in the Department of Radiology at the University of Chicago. His primary research interests are the objective assessment of image quality and image reconstruction in x-ray computed tomography.

Emil Y. Sidky received his BS degree in physics and mathematics from the University of Wisconsin, Madison, in 1987 and his PhD in physics from the University of Chicago in 1993. He held academic positions in physics at the University of Copenhagen and Kansas State University. He joined the University of Chicago in 2001, where he is currently a research associate professor. His current interests are in CT image reconstruction and large-scale optimization.

Xiaochuan Pan received his BS degree from Beijing University, Beijing, China, in 1982, and his MS degree from the Institute of Physics, Chinese Academy of Sciences, Beijing, in 1986, and his master’s and PhD degrees from the University of Chicago in 1988 and 1991, respectively, all in physics. He is currently a professor with the Department of Radiology, The University of Chicago. His research interests include physics, algorithms, and applications of tomographic imaging.

Citation Download Citation

Adrian A. Sánchez, Emil Y. Sidky, Xiaochuan Pan, "Use of the Hotelling observer to optimize image reconstruction in digital breast tomosynthesis," Journal of Medical Imaging 3(1), 011008 (22 December 2015). http://dx.doi.org/10.1117/1.JMI.3.1.011008

Submission: Received ; Accepted

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PDF + Save for Later

SHARE

GET CITATION

< Previous Article

|

Next Article >

KEYWORDS

Reconstruction algorithms

Sensors

Digital breast tomosynthesis

Signal to noise ratio

Image restoration

Signal detection

Optical filters

Create, Organize, Share

1.

Introduction

2.

Methods

2.1.

Digital Breast Tomosynthesis System Model

(1)

(2)

(3)

2.2.

Image Reconstruction Algorithms

(4)

(5)

2.2.1.

Filtered back-projection algorithm

(6)

(7)

(8)

2.2.2.

Back-projection filtration algorithm

(9)

(10)

2.2.3.

Λ-tomography

(11)

2.3.

Hotelling Observer Signal-to-Noise Ratio Metric

(12)

(13)

3.

Results

Table 1

Table 2

Fig. 1

Fig. 2

Fig. 3

Fig. 4

4.

Conclusions

Acknowledgments

References

Biography

Create, Organize, Share

Show All Keywords

KEYWORDS/PHRASES

SEARCH IN:

PUBLICATION YEARS

$Λ$ -tomography