2.1.

Mammograms

The mammograms we use were obtained from the UC Davis Medical Center and the University of Pennsylvania Medical Center under IRB-approved human-subject protocols that included deidentifying all image data. Images from 26 patients came from the UC Davis Medical Center and images from 9 patients came from the University of Pennsylvania Medical Center. Patient images were selected on the basis of having complete data, which consisted of two standard views (CC and MLO) of each breast for a total of 140 images. All of the mammograms were acquired on Hologic Dimensions systems (Hologic Inc., Bedford, Massachusetts) as part of breast cancer-screening exams. However, patient images in our sample were collected retrospectively on the basis of a recommendation for biopsy, and so they should not be considered a representative sample of the screening population.

A subset of 26 of these exams (acquired at UC Davis) was used previously in the original LFE paper.¹⁴ All breast image data were saved at the “for-processing” stage, after detector calibration had been performed, but before any display processing had been applied. This step was done to avoid any non-Gaussian statistics that might be introduced by transformations in the “for-display” processing.¹⁴ The images were log-converted from transmission to density before the analysis was performed. Pixel size in the mammograms was $70\mu \mathrm{m}$.

2.2.

Clustered-Blob Lumpy Backgrounds

The CBLBs have been developed by Bochud et al.^18,19 as an extension of Rolland’s and Barrett’s type I lumpy backgrounds²⁵ for applications in x-ray projection mammography. The idea behind lumpy backgrounds was that a background could be generated as a superposition of lump profiles placed at different locations throughout the background area. The CBLB extended this by allowing the lump profile to be a localized probability density for a cluster made up from smaller lumps.

A total of five previously published CBLB simulation procedures are considered here. These include the process described in the original publication¹⁸ and referred to as “OpEx99,” as well as four virtual breast phantoms developed in a subsequent publication.¹⁹ These are all sampled on a $100\text{-}\mu \mathrm{m}$-pixel size that was used in the publications. Two virtual phantoms based on oriented and isotropic processes are referred to as “SimpOri” and “SimpIso,” respectively, with “Simp” indicating that a single CBLB process is used. The remaining two phantoms both sum two CBLB processes modeling glandular and fibrous structures, respectively. These also incorporate oriented and isotropic processes and are referred to as “DoubOri” and “DoubIso,” with “Doub” indicating that a double CBLB process is used. The parameters for these processes are all taken directly from the original publications, with the only change being to scale the number of clusters to match the larger image size ($1024\times 2048$ instead of $256\times 256$) tested. For each phantom, 50 independent backgrounds are generated for analysis. Figure 2 gives examples of each of these processes [Figs. 2(a)–2(e)] along with a power-law Gaussian texture [Fig. 2(h)], sampled on a $100\text{-}\mu \mathrm{m}$-pixel grid, for reference.

Fig. 2

CBLBs and BT backgrounds. Samples of CBLBs are shown (a–e) along with those generated using BT processes (f, g). A Gaussian texture generated at $100\mu \mathrm{m}$ sampling (h) is also shown for reference.

It should be noted that the CBLB phantoms directly simulate a two-dimensional (2-D) image. This is considerably simpler than the methods below, which simulate a three-dimensional (3-D) breast and then project the 3-D simulation onto a 2-D detector.

2.3.

Binary Truncation Processes

The BT phantoms are based on the idea that the breast can be well modeled as mixture of adipose and glandular tissues, each characterized by its respective attenuation coefficient. These models of soft tissue contrast have typically neglected calcifications. The spatial distribution of the attenuation coefficients is generated by thresholding a 3-D random process, initially developed by Gong et al.²⁶ The approach of Abbey and Boone²⁰ (referred to as AB08) uses a 3-D Gaussian process truncated to $\sim 30\%$ volume glandular fraction. Reiser and Nishikawa¹⁰ (referred to as RN10) used phase randomization followed by thresholding to generate the distribution of adipose and glandular tissues and targeted a 75% glandular fraction.

It should be noted that both of these approaches have higher glandular fractions than the 19% average volume glandular fraction of the breast, as measured from breast CT images.²⁷ Both methods could be easily modified to meet this average glandular fraction by increasing the truncation threshold. There are other differences between the two processes (how they generate and truncate the 3-D random field) that can affect LFE. However, our interest is in evaluating whether LFE changes substantially between the different simulation procedures, and so we have adhered as closely as possible to the simulation procedures used in the publications.^10,20

Binary 3-D backgrounds were generated with an isotropic sampling of $100\mu \mathrm{m}$ in $x$, $y$, and $z$. The resulting grid was $1024\times 2048\times 512$, and thus approximately captured the 5-cm average thickness of a compressed breast in a screening mammography exam while achieving the size needed for a comparable region of interest (ROI) of the breast interior relative to FFDM images. Phantoms were projected onto the detector by summing the $z$-dimension of the 3-D phantom for a $1024\times 2048$ final image size. For each of the BT phantoms, a total of 50 independent backgrounds were generated for analysis. Examples of BT phantoms from these two processes are shown in Figs. 2(f) and 2(g).

2.4.

UPenn Virtual Breast Phantom

A virtual breast phantom development project has been underway at the University of Pennsylvania for several years.^21–23 This approach is based on simulating the major anatomical components of the breast, including skin, adipose tissue, fibroglandular tissue, and Cooper’s ligaments. The breast simulation is initially generated in an uncompressed state, to which a compression transform is applied using a finite element model.

We consider two versions of the UPenn phantom. One is the first published version of the phantom, tailored for projection mammography.²² In the second version, a new approach is used that involves simulating additional microstructure in the adipose compartments.²⁸ We refer to the original simulation procedure as the “coarse” phantom because of the strong texture of the projection images. The newer simulation is referred to as the “mixed” phantom because of the mixing of attenuation components in the adipose compartments. Mammograms are simulated from the 3-D phantom, assuming a 950-mL breast volume with 6.4 cm compressed breast thickness, a polyenergetic x-ray beam, and $70\text{-}\mu \mathrm{m}$ detector pixels similar to the Hologic detector. The phantoms are generated at an image size of $3328\times 4096\text{pixels}$. For each of the UPenn phantoms, 50 independent backgrounds are generated for analysis. Example images from these two phantoms are given in Fig. 3, along with a power-law Gaussian texture sampled on a $70\text{-}\mu \mathrm{m}$ pixel grid.

Fig. 3

Phantoms generated at $70\mu \mathrm{m}$ sampling. Samples from two versions of the UPenn virtual breast phantom are shown (a, b) along with a Gaussian texture (c) generated at $70\text{-}\mu \mathrm{m}$ pixel size and shown for reference.

2.5.

Image Power Spectra

This work makes extensive use of power spectra for both the sample images and the simulation procedures. Here, we briefly describe the windowed spectral estimation procedure we use. All image analysis (both power spectra and LFE) focuses on the ROIs within the image. Examples of the interior regions for the different phantoms are shown in Fig. 4. For the FFDM images, these ROIs are generated manually to indicate the interior of the breast, which is taken to be $\sim 10\mathrm{mm}$ inside the skin line of the mammogram and avoiding pectoral muscle, ribs, skin-folds, etc. The average area of the ROIs, across all patients and views, is $116.3{\mathrm{cm}}^2$ and ranges from 15.9 to $308.4{\mathrm{cm}}^2$. For the CBLB and BT phantoms generated on a $100\text{-}\mu \mathrm{m}$ grid and an overall simulation size of $1024\times 2048$, we used a central rectangle with an area of $112.3{\mathrm{cm}}^2$. This ROI is also used for the Gaussian process sampled on the same scale. For the UPenn phantoms (and corresponding $70\text{-}\mu \mathrm{m}$ Gaussian phantom), the ROI is also generated similar to the FFDM images. Since the UPenn phantoms all have the same breast contour with randomized interior structure, a single ROI ($106{\mathrm{cm}}^2$) is appropriate for all of these phantoms.

Fig. 4

ROIs used for analysis. Each image shows a phantom with the ROI indicated. (a) The CBLBs and BT processes used a rectangular ROI with an area of $112{\mathrm{cm}}^2$. (b) The UPenn phantoms and (c) FFDM patient images used manually determined the ROIs. The UPenn phantom ROI has an area of $106{\mathrm{cm}}^2$. The ROI in the FFDM image shown (LCC view) is roughly the average size of the area in the image set ($111{\mathrm{cm}}^2$). In this image, a small skin fold in the lower right corner is explicitly avoided in the ROI.

The power spectrum is estimated from square image patches that are randomly selected so that the entire patch falls within the interior of the breast. Let $g_i[n,m]$, $n,m=0,\dots ,N-1$, be the sample patch from the $i$’th FFDM image or the $i$’th image generated by one of the simulation procedures ($i=1,\dots ,I_{\mathrm{Samp}}$). The patch size is fixed at $512\times 512\text{pixels}$ (i.e., $N=512$). For the FFDM images where the breast was too small to fit, the patches are excluded from the power-spectrum estimate. This constraint leaves 34 cases for use in each of the MLO views, and 33 cases for each of the CC views. The simulations all use 50 sample images. A windowed power spectrum is estimated for each ROI in a sample image as

We use radial averages to convert the 2-D power spectra in Eq. (1) to a one-dimensional function of radial frequency, consistent with previous works.^4,6,30 This has typically been justified on the basis of approximate isotropy of power spectra for breast images, although other investigators have shown orientation dependence.³¹ Since our study is focused on the power spectrum at a more general level, we use radial averages and neglect any orientation dependence in our data. Let $f_{k,l}$ represent the radial frequency for the indices, then

Let $q$ be the radial frequency index ($q=0,\dots ,N/2$) with associated radial frequency ${\rho }_q=q/N\mathrm{\Delta }$. We define $F_q$ as the set of all 2-D frequency index pairs in a band centered on ${\rho }_q$ that has a width of one radial frequency sample (i.e., $F_q=\{[k,l]:-1/2N\mathrm{\Delta }\le f_{k,l}-{\rho }_q<1/2N\mathrm{\Delta }\}$); let $N_q$ be the number of index pairs in the band, then the radial power spectrum is given as

The radial power spectrum for the process is generated by averaging over the image samples. In addition, we find it useful to work with the log-power spectrum at times. We generate the average log-spectrum by averaging over samples as well

The sample standard deviation across $i$ is used to quantify the variability of the log-power spectrum

2.6.

Simulating Acquisition Noise

The above simulation procedures describe how the anatomical components of the phantom images are generated. However, a full simulation should also include the generation of acquisition noise. In the initial conference report of this study,²⁴ acquisition noise was modeled as white Gaussian noise. While this had the correct qualitative effect of flattening the power spectrum at high frequency, it is noticeably different than the actual power spectra of digital images. In addition, the different simulation procedures are normalized differently leading to large differences in the amplitude of the power spectra.

In this study we model the acquisition noise as another power-law process that is fitted to the power spectrum of the FFDM images over a frequency range from $0.1\mathrm{cyc}/\mathrm{mm}$ to the Nyquist limit. The noise process has two free parameters: a power-law amplitude, $A_{\mathrm{Nse}}$, and exponent, ${\beta }_{\mathrm{Nse}}$. We also scale the simulation processes by a factor $A_{\mathrm{Scl}}$ to have approximately the same amplitude as the FFDM images. These effects are achieved by fitting the power spectra of the simulation procedures to the power spectrum of the FFDM images. Let ${\widehat{S}}_{\mathrm{Bck}}[k,l]$ be the power spectrum estimated for one of the background simulation procedures. We fit the simulation power spectrum to the image spectrum using a three-parameter model

2.7.

Laplacian Fractional Entropy

LFE is based on the relative entropy ($K\text{-}L$ divergence) between the response histogram of a Gabor filter, normalized to be a probability distribution, and a best-fit Gaussian distribution over the same histogram bins. If the histogram represents an underlying Gaussian process, then it will be well fit by the Gaussian distribution, and the relative entropy will be negligible. Any sort of departure from a Gaussian distribution in the histogram will lead to nontrivial relative entropy. This relative entropy is normalized by the relative entropy of the best-fit Laplacian distribution with respect to this Gaussian. In this way, the Laplacian relative entropy serves as a yardstick for interpreting the Gaussian relative entropy. As a result, an LFE value of 100% means that the histogram exhibits as much departure from a Gaussian profile as the Laplacian distribution does. In principle, any non-Gaussian distribution could be used in place of the Laplacian. However, Laplacian distributions have played an important role in the analysis of natural scenes.^32–34

The LFE was computed according to the procedure specified in a previous publication.¹⁴ Gabor filters that spanned 11 center frequencies from 0.125 to $4.0\mathrm{cyc}/\mathrm{mm}$ were evaluated at six different orientations (0 deg, 30 deg,…, 150 deg). The filters were sine-phase with a bandwidth of 1.4 octaves and an aspect ratio of 1. Filter responses were generated by convolving each filter with an image. All responses from 1 cm inside the boundary of the ROI were used to form response histograms from which the LFE measures at different frequencies were estimated. The histograms binned the central 99% of the responses, with an additional 1% bin for the remaining extremal values. The computational scheme for determining LFE for a particular Gabor filter in a given image is shown in Fig. 5, resulting in an LFE value of 63%.

Fig. 5

Computation of LFE. The image (a) is convolved with a Gabor filter (upper left corner), and all filter responses in the ROI are retained (b). (c) The resulting histogram is normalized to a probability distribution by dividing by the total number of responses and is shown with fitted Gaussian and Laplacian distributions. The LFE measure is the relative entropy of the data to the Gaussian divided by the relative entropy of the Laplacian to the Gaussian, which is 63% in this case.

For the FFDM images, LFE estimates from a single filter are averaged across the six orientations and across the four views for each patient. The result is averaged over patients, to get a final estimate of LFE at each spatial frequency. The sample standard deviation over patients is used to quantify the variability of LFE. For the simulation processes, we average across orientations and samples, to get an estimate of mean LFE for the process. The standard error of the mean (across samples) is typically quite small ($<4\%$ of the mean on average).

3.1.

Power Spectra

A comparison of radially averaged power spectra for simulation procedures and FFDM images are shown in Fig. 6. The simulation spectra have been fitted using the procedure described in Sec. 2.6. Each plot shows the results for a different class of simulation procedures. The fitted amplitudes of the acquisition power spectra ranged from 156.8 to 371.5, while the power-law exponents of the acquisition noise models ranged between $-0.27$ and $-0.90$.

Fig. 6

Radial power spectra. Radial averages of the power spectrum are shown for each of the simulation approaches. Each of the phantom procedures has been scaled and has Gaussian noise added to match the observed FFDM power spectrum across a frequency range of 0.1 to $5\mathrm{cyc}/\mathrm{mm}$. The FFDM data are shown on each plot for reference with error bars of $\pm 1\mathrm{sd}$ for comparison.

In most cases the fitting procedure results in good agreement between the fitted model and the FFDM data. The average relative error across all frequencies in all the simulation procedures is 6%. The largest divergence between fitted model and data occurs for the OpEx99 model of the CBLB [Fig. 6(a)] at low spatial frequencies ($<0.3\mathrm{cyc}/\mathrm{mm}$) where relative error is over 100%, but these drop below 10% at higher frequencies. Fitted spectra from the binary processes [Fig. 6(b)] and the UPenn breast phantom [Fig. 6(c)] agree with the sample estimates across the tested frequency range. The Gaussian processes used for reference [Figs. 6(b) and 6(c)] also agree well across the spectrum.

Measured from spatial frequencies between 0.1 and $1.0\mathrm{cyc}/\mathrm{mm}$, the power-law exponent of the FFDM power spectrum is $\beta =-3.01$, which is consistent with previously published results.^6–8 With the exception of the OpEx99 process ($\beta =-3.81$), the relative errors in the power-law exponents are within 10% of the sample estimate. The Gaussian processes, generated with a power-law slope of 3.0, result in relative errors of $<3\%$, which suggests that the spectral estimation procedure is not introducing a substantial bias in the power-law exponent.

Even with some departures at low frequencies by the CBLBs, the general finding is that the power spectra of the virtual phantoms are reasonably well fit by a power-law over the frequencies in which anatomical variability dominates, even though the various phantoms have quite different visual appearances. This is a testament to the ubiquitous nature of the power-law spectrum. However, when the goal is to distinguish between various models, the generality of the power law may work against its usefulness. This serves as the motivation for investigating higher-order statistical properties.

3.2.

Laplacian Fractional Entropy

LFE results for the various simulation procedures are shown in Fig. 7. Here again, the FFDM results for left CC view images are plotted (with $\pm 1\text{sd error}$ bars across patients) for comparison with each class of simulation procedure. Other views have similar LFE values. The general profile of the LFE curve is similar to previous reports,¹⁴ which is not surprising, given that 26 of the 35 cases were used for those results. LFE appears to be dropping across the lower frequencies with a mild shoulder around $1\mathrm{cyc}/\mathrm{mm}$, where acquisition noise begins to dominate, driving the measure to zero.

Fig. 7

LFE plots. Plots are shown for (a) the CBLBs, (b) BT backgrounds, and (c) two versions of the UPenn virtual breast phantom. The LFE data for the patient FFDM images are shown on each plot for reference, with error bars of $\pm 1\mathrm{sd}$ for comparison. In the binary processes plot (b), there is also a plot from the Gaussian process generated on a $100\text{-}\mu \mathrm{m}$ pixel size for reference and, on the UPenn virtual breast phantom plot, there is also a plot for the Gaussian process generated on a $70\text{-}\mu \mathrm{m}$ pixel size.

Figure 7(a) shows LFE plots for the CBLBs, along with a matched Gaussian texture. These plots all have fairly small values across the frequency range ($<5\%$). These are similar to the matched Gaussian process plotted in Fig. 7(b). Nonzero values of LFE for the Gaussian textures reflect departures from the Gaussian profile in the histogram due to the limited spatial extent of the ROI used to compute the histogram and the finite sample of histogram responses.

None of the CBLBs show any substantial non-Gaussian structure by the LFE measure. This is somewhat surprising given that the procedure for generating the phantoms is rigorously non-Gaussian.^18,19 However, it appears that—seen through the lens of Gabor receptive fields—there is little higher-order structure in these backgrounds. Given that the number of background “blobs” in a simulation are generated using sums over Poisson random deviates, it may be that for the frequencies we are testing, Poisson variability is well approximated as Gaussian. It is also possible that substantial non-Gaussian structure in the images is not being captured in the Gabor-filter histograms used by the LFE measure. Further analysis is needed to resolve this issue.

Figure 7(b) shows LFE plots for the phantoms based on the BT processes, along with the matched ($100\mu \mathrm{m}$) Gaussian texture and the FFDM samples. The AB08 process results in higher LFE than the RN10 process, and this is mostly due to more regions of uniform intensity, which add to the response histograms near zero. This explains the increased LFE observed here compared to the earlier conference report.²⁴ Future implementations and optimizations of these processes should carefully consider the simulation parameters including sampling resolution and approaches to dealing with power-law processes near $f=0$.

Figure 7(c) shows LFE plots for the UPenn phantoms along with the matched ($70\mu \mathrm{m}$) Gaussian texture, and the FFDM sample estimates. As expected, the Gaussian texture has very low LFE values similar to the $100\mu \mathrm{m}$ Gaussian texture plotted in Fig. 7(b). This serves as a check that the different pixel sizes do not produce any unexpected results. The coarse UPenn phantom shows a substantial peak in LFE at frequencies between 0.3 and $1.5\mathrm{cyc}/\mathrm{mm}$ that is not present in the FFDM data. This peak is substantially modulated in the mixed UPenn phantom leading to a better agreement with the FFDM data over the entire spectrum.

The UPenn phantom results show how LFE can be used to guide the development of virtual breast phantom procedures. Both the coarse and the mixed phantoms can be made to closely match the power spectrum of breast images, which means that the power spectrum is not a useful way to identify which one is closer to the real images. In contrast, LFE shows a substantial difference between the two and suggests that mixing in the adipose compartments is moving the phantom in a more realistic direction. This is not entirely surprising given that there are physiological reasons to expect a distribution of attenuation values within a glandular compartment, but it is advantageous to have a quantitative measure that reflects this.

As a way to quantify how closely the phantom LFE matched the FFDM LFE at frequencies of interest, we compare differences in the average phantom LFE and FFDM LFE relative to the standard deviation in LFE across patients in the FFDM data. For frequencies between 0.125 and $1\mathrm{cyc}/\mathrm{mm}$, the average error of the AB-08 phantom relative to the standard deviation of the FFDM LFE is 52%, and for the mixed UPenn phantom it is 125%. For the RN-10 phantom, the coarse UPenn phantom, all the CBLB phantoms, and the Gaussian processes, the errors are $>200\%$.

The LFE results presented here are somewhat different than in our initial results.²⁴ This demonstrates the potential impact of varying the simulation procedures. For this reason, and others noted below, we hesitate to accept this as a reflection of the general validity of the virtual phantom procedures.

3.3.

Study Limitations

Our study is fairly limited, and hence we caution against the overinterpretation of our results, and in particular the establishment of a “best” phantom on the basis of these results. The purpose of our study was to demonstrate that phantoms which have similar power-law spectra may exhibit considerable differences in higher-order statistics, and that LFE can be used to quantify those differences in a meaningful way. Comparison with the LFE derived from FFDM images from a limited set of patients scheduled for biopsy was made to give a general sense of how real FFDM images behave. This study was not designed to identify the most appropriate phantom. A more focused evaluation of phantom approaches should consider the population of patients from which sample images are used, and it would likely also be useful to allow for some training of parameters in the phantom procedure as well.