In this section, we first present several phantom data reconstructions to illustrate the robustness of the imaging system and the algorithm. We used a homogeneous phantom to test for image artifacts that might result from systematic errors, calibration, geometric modeling, and algorithmic inaccuracy. The inhomogeneous phantoms, on the other hand, were used to test the sensitivity of the imaging data to objects that mimic the size and contrast of a tumor. A range of patient data reconstruction results follow and are shown in four groups based on the sizes of the glandular regions. For each group, the common findings are summarized and then further compared across groups. The mean bulk optical properties of each group were tabulated to facilitate comparisons to other published studies [42
For all patient reconstructions, a finite element forward mesh with uniform node spacing of 4 mm and a reconstruction mesh with node spacing of 8 mm were generated from DBT images. This typically results in 10 000-50 000 forward nodes and 1500-6000 reconstruction nodes depending on the volume of the breast. The FEM forward solver uses an optimal block size of 8 [18
] to solve the solutions for all sources and detectors by an iterative block solver [19
]. On a Pentium Xeon 2.4-GHz CPU PC, the average solving time for a solution with 20 000 unknowns is about 1 s.
The source/detector positions for phantoms and clinical data reconstructions are plotted in . The source/detector coverage for clinical measurements were specially optimized using a histogram generated from DBT images from 40 patients to capture 80% of the breast coverage. For phantom reconstructions, only RF data were included, therefore, we only plot the RF source/detector positions in . In all cases, both the source and detector probes are parallel to each other and directly contact the target surface through a thin transparent acrylic paddle.
Illustrations of sources/detector configurations: (a) sources and (b) detectors for phantom measurements; (c) sources and (d) detectors for patient measurements.
The Gauss-Newton iterations were applied to both the bulk property estimation and the full-image reconstruction, where the maximum iteration numbers were set to 20 and 5 for phantom and clinical data, respectively. In both cases, the SD coefficients and optical properties were estimated simultaneously. As we discussed in Section II-D5, we used the estimated bulk optical properties and the averages of the estimated SD to set the homogeneous initial guess for the full-image reconstruction. Both the RF and CW MUX measurements were used in the bulk estimation and image reconstruction of the patient measurements to achieve improved image resolution and accuracy. For a typical reconstruction, constructing the Jacobian matrix takes about 2 s, with another 4-5 s needed to solve the update equation. The average iteration time for the reconstruction ranges from about 1-3 min (depending on the number of useful sources and detectors).
A. Phantom Reconstructions
The phantom used in this paper is a balloon filled with intralipid/ink solution. The balloon is vertically compressed between the two compression paddles to mimic a breast under compression. The dimensions of the compressed balloon are roughly 15 cm × 18 cm horizontally and 5 cm in the compression direction. A picture of the compressed balloon is shown in . The liquid solution inside the balloon contains 0.8% intralipid and 0.017% Higgins India ink (Stanford, IL). The horizontal extent of the balloon was identified by tracking the contours from the top-view photo in and mapping to the probe coordinates assisted by landmarks. Smooth nonuniform rational B-spline (NURBS) patches were used to connect these contours, from which the surface and volumetric forward/reconstruction meshes were generated .
Balloon phantom experiment illustration: (a) top-view photo of the compressed balloon, (b) forward FEM mesh, and (c) inclusion.
The raw measurement of the balloon phantom was calibrated by (16
). In this example, we used the 830-nm RF measurement only. The bulk properties of the balloon phantom estimated from (23
) are μa
= 0.063 cm-1
at 830 nm. An independent measurement on a similar sample with a 8-wavelength RF NIR spectrometer (ISS Inc., Champaign, IL) was performed and the measured optical properties are μa
= 0.081 cm-1
at 830 nm. Given the significant differences between the two instruments, including measurement geometry, probe design, wavelengths, the discrepancy between the two estimates seem comparable with published studies on cross-instrument validation [50
]. Full-image reconstructions of the balloon phantom yield the absorption and scattering profiles at 830 nm. Horizontal slices extracted from 3 cm above the bottom of the phantom are shown in . The reconstructed images demonstrate reasonable homogeneity in the target domain with maximum perturbations of ±2.1% in the absorption and ±1.9% in the scattering image.
Fig. 7 Reconstructed image slices for a balloon phantom filled with intralipid solution: (a)μa (cm-1) image and (b) image. Both images were extracted at the horizontal plane 3 cm from the bottom of the phantom. Maximum perturbations for the absorption (more ...)
A 2.5-cm-diameter glass sphere filled with the same solution except with doubled ink concentration was positioned inside the balloon and manipulated by moving a wood-stick through a translation stage. A picture of the inclusion is shown in . We first position the sphere near the center of the balloon, and then moved it slightly to the right side. Vertically, the sphere is close to the bottom plate. In both cases, top-view pictures were taken while vertically lifting the inclusion in order to track its position. With the measurement calibrated by a homogeneous balloon measurement, we performed optical image reconstruction for the two measurements and slices of the absorption images are shown in . A 2.5-cm-diameter circle was overlayed on top of the images based on the position from the top-view photos in both cases.
Fig. 8 Reconstructed μa (cm-1 images of a balloon with a 2.5 cm diameter spherical inclusion: (a) configuration 1 and (b) configuration 2. Image slices were extracted from a horizontal plane 1.5 cm above the bottom plane. Contour of the inclusion is (more ...)
From these images, we can see that the positions and the sizes of the target were recovered fairly accurately. However, the maximum absorption contrast in both cases are 1.3:1 and 1.5:1, respectively, which are both lower than the expected contrast of 2:1. The underestimation of the contrast was considered a result of the smoothing effect of the regularization and as previously reported in similar reconstructions [52
]. The thin glass wall of the inclusion may also contribute to the underestimation; however, given the contrast reduction due to smoothing of the image reconstruction as shown for example in [52
], the systematic error introduced by the glass-wall may be negligible.
B. Clinical Data Reconstructions
With confidence in both the measurement data and reconstruction algorithm, we applied the above analysis to the clinical data acquired from healthy subjects between August 2005 and December 2006. The processing steps for the clinical data are generally identical to those for the phantoms, with the exception that the forward and reconstruction meshes were generated from 3-D tomosynthesis images with the registration and mesh generation described in Section II-D2.
A total of 49 healthy subjects were enrolled in our study, including cases of single-side and bilateral breast measurements for a total of 68 data sets with both valid optical and DBT measurements. A semi-automatic segmentation software, Insight-SNAP [53
], was used to generate the 3-D volumetric segmentation of the fibroglandular regions from the DBT scans, from which a fibroglandular volume fraction (g
) was calculated by taking the ratio of the fibroglandular volume by the imaged breast volume. Using this metric, we divided all imaged breasts into groups and the sample size of each group, the means and standard deviations of the recovered bulk oxygenated-hemoglobin concentration (HbO), deoxygenated-hemoglobin concentration (HbR), total hemoglobin concentration (HbT = HbO + HbR), oxygen saturation (SO2
= HbO/HbT), and reduced scattering coefficient
at 830 nm are listed in . From this table, we can see that with the increase in g
, an increase in HbT, HbO, and HbR is also observed except for the last group, while the SO2
values show a decreasing trend. The mean values of HbT is relatively low compared with literature values for uncompressed breasts [42
]. However, due to the applied mammographic compression, removal of blood is anticipated from our previous studies [25
] and the values in are reasonable.
TABLE I Average Bulk Properties for 68 Healthy Breasts. g IS the Fibroglandular Volume Fraction; N is the Number of Breasts for Each Group. The Units are Micrometers for HbO, HbR, and HbT and cm-1 for
Twelve DBT image slices (six for LMLO view, six for RMLO view) are shown in (right-hand side image of each image pair). In these images, one can easily identify the boundaries of the breast, the fibroglandular tissue (circled by thin contour lines), the fatty tissue (areas surrounding the fibroglandular region), and chest-wall muscle [the triangular regions delimited by white dashed lines in cases and ].
Fig. 9 Reconstructed HbT and DBT images for various subjects. The images in the first row [(a)-(c)] are examples of fatty breast reconstructions; the images (d)-(f) scattered density breasts; images (g)-(i), heterogeneously dense breasts; and (j)-(l), extremely (more ...)
Reconstructed HbT and DBT images for 12 representative healthy breasts are shown in . They are divided into four groups based on their radiological breast densities: fatty (FA), scattered density (SC), heterogeneously dense (HD), and extremely dense (ED) [54
]. For each group, we show glandular region outline overlays on top of the HbT and X-ray images. A thick dashed line was used to mark the field-of-view of the optical probes. The optical images between the dashed line and the nipple are not expected to have reconstructed contrast because of the extinction of measurement sensitivity. All slices were extracted from the horizontal plane 2 cm above the bottom of the compressed breast. The HbT values for the images were computed from absorption coefficients at 685 and 830 nm by fitting the oxygenated/deoxygenated hemoglobin absorption spectra to the data while assuming water and lipid concentrations are 23% and 58%, respectively. The water and lipid concentrations were reasonable extrapolations from literature values 31% and 57%, respectively [42
], considering a reduction due to mammographic compression. To investigate the sensibility to water and lipid, we doubled water concentration and found that the HbO values decreased only about 20% and the HbR values increased less than 10%. Reducing the lipid concentration to 20% made less than 2% changes to HbO and HbR.
Here, we summarize the observations by examining the images for each group in . For breasts with very little glandular content, i.e., , the averaged HbT values are relatively small compared with the other three groups, similarly for the perturbations in HbT distributions. For breasts with a small amount of glandular tissue, i.e., those in , a positive contrast in HbT from the bulk values at the glandular region appears in the images. The low-HbT region is located mostly on top of the fatty tissue. The triangular chest-wall muscle areas for all three cases were successfully reconstructed and present the highest HbT concentrations in the image. For heterogeneously dense breasts, the findings are largely identical to the scattered density group, except for , where the positive contrast of the glandular regions appears slightly deformed by the low-HbT region at the top-right of the breast. For breasts with a high percentage of glandular tissue, the image features largely resemble those of the fatty breast group, i.e., relatively small contrast in the breast; however, the baseline HbT for this group is about twice as high. We also noticed that the contours of the recovered glandular regions in the dense breast group contain a somewhat low HbT region similar to . The HbT decrease in these regions ranges from about 5% to about 20%.
To generalize the findings from , we performed a region-of-interest (ROI) analysis for the breasts with small to moderate amounts of fibroglandular tissue and simultaneous RF and CW measurements (N = 26). The ROIs were manually created based on the X-ray segmentation and the field-of-views of the optical probes for each breast. The average HbT of the fibroglandular tissue is 20.1 ± 6.1 μm (N = 24), that for adipose is 15.4 ± 5.0 μm (N = 26), and for muscle 22.2 ± 7.3 μm (N = 17). Paired t-test indicates that the HbT of the fibroglandular tissue in these breasts is significantly greater than that of adipose tissue (N = 24, p < 0.0001) with an average contrast of 1.4:1. Note that the estimated contrast is likely lower than the true contrast as shown in our phantom results.
distributions are similar to those of HbT for each patient. Representative images of the
from two breasts, i.e., , are shown in , respectively, where the glandular regions demonstrate higher scattering values. The corresponding SO2
images of the two subjects  seem to present lower values within the fibroglandular segmentations, which suggest higher metabolic activities in the fibroglandular tissue.
Fig. 10 Image slices for (a) at 830 nm, (b) oxygen-saturation (SO2) corresponding to case (i) in Fig. 10, and (c) at 830 nm, (d) SO2 corresponding to case (d) in Fig. 10. The markings on the images are defined similarly as in .
We can conclude from these observations that: 1) the chest-wall muscles, where visible, normally present higher HbT and sharper edges in the image; 2) the glandular regions within the coverage of the optical detectors have a roughly 10% increase in the HbT images; 3) both fatty and dense breasts show less variation within the image than the heterogeneous group; 4) for some cases, a region of slightly decreased HbT can be seen partially overlapping on the glandular area around the central part of the breast.
C. Pressure Simulations
The presence of a decreased-HbT region in some images is a unique feature for this study and deserves further explanation. As we know, mammography scans are typically associated with applying large compression to the breast. It is reasonable to assume that the blood in the tissue is redistributed, driven by the gradient of tissue internal stress. In this case, we would expect lower-HbT content in the region where tissue stress is stronger.
The key component to validating the hypothesis above is to obtain the stress distribution and correlate it with optical images. To accurately model the breast under large deformation is fairly complicated and beyond the scope of this paper. In this section, a mechanical simulation is presented under simplified assumptions to qualitatively support our hypothesis.
The simulation was performed with multiphysics software ANSYS ED. The simulated breast geometry is shown in . It was created by cutting a sphere with a radius of 8 cm centered at the origin by three planes at z
= -3 cm, z
= 3 cm, and y
= 2 cm, respectively, and subsequently scaling the model in the z
axis by a factor of 5/6. The breast model was treated as a linear isotropic material with Young's modulus of 25 MPa and a Poisson constant of 0.48 (for semi-volume-preserving) [56
]. Two displacement boundary conditions were applied to the upper and bottom surfaces, respectively, with Δz
= -0.2 cm at the top and Δz
= 0.2 cm at the bottom. No constraints were applied for the x
directions for the two flat surfaces and other surfaces.
Breast compression simulation and von Mise's stress crosscuts: (a) breast model, (b) (x, y) cross section, (c) (y, z) cross section. The darker the color scale, the higher the stress.
The von Mises stress plots on the (x,y) and (y,z) cross sections of the compressed breast are shown in . As hypothesized, the center part of the breast has higher stress than the periphery regions. Driven by the gradient of the internal stress, the blood in breast tissue will be “squeezed” away from the high-stress regions and results in the low HbT regions in the images shown in the previous section.
One must note that the simulation shown above presents a simplified view of the effect of compression on breast tissue. These static features of stress distribution provide a qualitative justification for the patterns of hemoglobin distribution observed in our reconstructed images. However, the establishment of pressure gradients due to external compression sets in motion several other processes [25
Immediately after breast compression, external force leads to an instant increase in interstitial fluid pressure in the breast sufficient to squeeze some of the blood out of the tissue and consequently leads to an immediate drop in overall blood volume. At the same time, breast tissue (and biological tissue in general) is visco-elastic [56
], meaning that the stress-strain response curve of breast tissue is different from that of linear elastic materials in that it exhibits “creeping” (relaxation) due to water movement and collagen matrix reorganization [57
]. Also, the breast has a heterogeneous internal structure composed of both higher stiffness fibroglandular tissue as well as less stiff adipose tissue. Furthermore, breast tumors have been shown to have even higher stiffness [59
] and interstitial pressure [60
] than normal breast tissue. The combined effect of these features is a dynamic stress redistribution following compression manifested by a partial return of blood to tissue and stress/strain transfer from the stiffer areas to the less stiff ones.
The additional complexities of the tissue compression response also represent an opportunity for identifying dynamic characteristics useful in cancer diagnosis. However, this is beyond the scope of the current publication.