Combined Optical Imaging and Mammography of the Healthy Breast: Optical Contrast Derived From Breast Structure and Compression (original) (raw)

• Journal List
• HHS Author Manuscripts
• PMC2642986

IEEE Trans Med Imaging. Author manuscript; available in PMC 2009 Feb 15.

Published in final edited form as:

PMCID: PMC2642986

NIHMSID: NIHMS85682

G. Boverman is with the is with the Department of Biomedical Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180 USA (e-mail: ude.ipr@grevob).

D. B. Kopans and R. H. Moore are with the Avon Foundation Comprehensive Breast Evaluation Center, Massachusetts General Hospital, Boston, MA 02114 USA (e-mail: gro.srantrap@snapokd; gro.srantrap@eroomr).

E. L. Miller is with the Department of Electrical and Computer Engineering, Tufts University, Boston, MA 02155 USA (e-mail: ude.stfut.ece@rellimle).

D. H. Brooks is with the Department of Electrical and Computer Engineering, Northeastern University, Boston, MA 02115 USA (e-mail: ude.uen.ece@skoorb).

*Q. Fang is with the Massachusetts General Hospital, Charlestown, MA 02148 USA (e-mail: ude.dravrah.hgm.rmn@qgnaf).

Abstract

In this paper, we report new progress in developing the instrument and software platform of a combined X-ray mammography/diffuse optical breast imaging system. Particularly, we focus on system validation using a series of balloon phantom experiments and the optical image analysis of 49 healthy patients. Using the finite-element method for forward modeling and a regularized Gauss-Newton method for parameter reconstruction, we recovered the inclusions inside the phantom and the hemoglobin images of the human breasts. An enhanced coupling coefficient estimation scheme was also incorporated to improve the accuracy and robustness of the reconstructions. The recovered average total hemoglobin concentration (HbT) and oxygen saturation (SO2) from 68 breast measurements are 16.2 μm and 71%, respectively, where the HbT presents a linear trend with breast density. The low HbT value compared to literature is likely due to the associated mammographic compression. From the spatially co-registered optical/X-ray images, we can identify the chest-wall muscle, fatty tissue, and fibroglandular regions with an average HbT of 20.1±6.1 μm for fibroglandular tissue, 15.4±5.0 μm for adipose, and 22.2±7.3 μm for muscle tissue. The differences between fibroglandular tissue and the corresponding adipose tissue are significant (p < 0.0001). At the same time, we recognize that the optical images are influenced, to a certain extent, by mammographical compression. The optical images from a subset of patients show composite features from both tissue structure and pressure distribution. We present mechanical simulations which further confirm this hypothesis.

Keywords: Breast imaging, multimodality imaging, tomography

I. INTRODUCTION

BREAST cancer is among the top threats to women's health. In the U.S., breast cancer accounts for 15% of cancer mortality and ranks second in all cancer deaths in the female population [1]. Early detection of breast cancer is critical to reduce the mortality rate [2]. Several novel noninvasive imaging methods have been proposed and applied to breast cancer early diagnosis. Near-infrared (NIR) spectroscopy/tomography [3]-[5] shows great promise for breast cancer detection and treatment monitoring [6]-[13]. In this imaging modality, NIR light is strongly scattered by tissue structures and absorbed by chromophores as it propagates through the breast tissue. In the 650-900 nm spectral range, the important chromophores are oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR). By performing optical measurements at multiple wavelengths, the HbO and HbR concentrations can be extracted by fitting measured light attenuations using the known tissue chromophore absorption spectra. In diffuse optical tomography (DOT), the absorption spectra fitting is performed not only on local measurements, but also on the measurements between spatially distributed source and detector arrays. In this case, 2-D or 3-D images of HbO and HbR concentration maps can be produced [7], [10], [14]. Because the tumors normally present metabolism and vascular structure significantly different from those of normal tissues, the hemoglobin images from NIR imaging are potentially useful in assessing tissue status and malignancy [15], [16].

Because of the nonlinear nature of the light diffusion and subsequent severe illposedness of the inverse problem [14], diffuse optical tomography generally produces smooth images with spatial resolution much lower than that of the standard X-ray mammography or the emerging tomosynthesis (DBT) technique. However, the ability to image tissue metabolism using NIR measurements is rather attractive for assisting diagnosis. To this end, we have implemented a multimodality approach: combining NIR tomography with X-ray mammography to produce accurately co-registered optical and X-ray images. Because of the spatial co-registration, the optical images, i.e., HbO, HbR, and oxygen saturation (SO2), can be better interpreted in conjunction with high resolution X-ray images; on the other hand, the inherent drawbacks of mammography, such as low contrast between lesion and normal tissues, especially in dense breasts, and susceptibility to unclear lesion boundaries, can be partially overcome with the additional functional information from NIR imaging, resulting in improved diagnosis. Functional overlay on structural images in this combined approach is analogous to that in the combination of positron electron tomography and computed tomography (PET/CT) [17].

A combined X-ray/optical breast imaging system was built at Massachusetts General Hospital (MGH), Charlestown, between 2001 and 2005. The original prototype of the system included only a radio-frequency (RF) modulated imaging unit [10]. After the first clinical trial of the system in 2004, significant modifications were applied. These changes include an additional Galvo multiplexer (MUX) with three continuous-wave (CW) lasers, a dedicated CW unit for continuous monitoring of tissue status at two wavelengths, an optimized RF laser wavelength spacing, and more flexible probe design. The additional multiplexed CW lasers greatly improved the measurement spatial sampling, and potentially improved the image resolution; the new 685-nm RF lasers and the extended CW wavelengths (Fig. 1) also helped to determine the hemoglobin concentration more accurately. Note that the selected wavelengths do not have sufficient sensitivity to detect tissue water and lipid concentrations, therefore, the recovery of the hemoglobin distribution is the major focus of this paper.

Block diagram of TOBI breast imaging system. The physical distributions of source and detector fiber positions can be found in Fig. 5.

In parallel, the optical image reconstruction algorithms have been improved dramatically in both accuracy and computational efficiency. A finite element method-based forward solver [18] was developed to replace the Born approximation with infinite-slab assumption used previously [10]. The use of an FEM solver allows the utilization of an unstructural mesh to model the curved breast boundary extracted from DBT images. It is also powered by a highly optimized iterative multiright-hand-side solver and is able to compute multiple forward solutions with a single forward-backward substitution [19]. A Gauss-Newton iterative reconstruction algorithm [14], [20] was implemented to produce reliable images for measurement data even when corrupted by large source/detector coupling variabilities (or SD) [21]. A detailed discussion on the forward and inverse algorithms can be found in Section II.

A clinical trial of the combined X-ray/optical imaging system started in August 2005. To date, we have scanned 65 subjects with the improved system and software platform (39 of them have bilateral breast measurements). Forty-two of those subjects were recruited from screening patients and 23 from patients recalled for diagnosis. The diagnostic results indicated that 49 patients were healthy, 12 had tumors, and four had benign lesions. In this paper, we focus on image reconstruction and analysis of the data from healthy patients. The primary goals of this analysis are 1) to outline the key components and improvement in a working prototype of a combined X-ray and DOT imaging system, 2) to characterize the system performance by phantom experiments, 3) to interpret the reconstructed breast images by directly comparing the structural and functional images produced by the system, and 4) to use the findings to guide the future development of multimodal DOT breast imaging.

In Sections II and III, we discuss the experimental procedures and analysis workflow for image reconstructions of healthy subjects using our imaging hardware and software platform. Sections II-A and II-B focus on the imaging system and the experiment protocols. They are followed by five more subsections expanding on the algorithmic details, including the forward/inverse models and coupling coefficient estimation. In Section III, we first report phantom imaging examples to assess the performance of the system. For clinical data reconstructions, we tabulate the average bulk optical properties for all subjects based on the density of their breasts. We also report the reconstructed images of 12 representative healthy breasts with region-of-interest overlays, and interpret the images by comparing with the corresponding X-ray scans. Our conclusions are summarized in Section IV, together with a discussion of the prospects for this combined X-ray/DOT approach.

II. METHODS

This section is dedicated to elaborating the experimental settings and image reconstruction details for this analysis. The first half of this section, including Sections II-A-II-C, focuses on the clinical aspect of the study. Readers can find an overview of the combined imaging system, the experimental protocols and data analysis steps in this section. The second half of this section, including all subsections of Section II-D, further discusses forward modeling, image reconstruction algorithm and measurement calibration.

A. System Overview

A block diagram of the tomographic optical breast imaging system (TOBI) is shown in Fig. 1. A picture of the combined TOBI/tomosynthesis (DBT) system in the clinical environment can be found in Fig. 2, shown with a close-up view of the optical probe (Fig. 3).

Picture of the combined TOBI/DBT system in clinical environment. The TOBI system include both RF and CW source/detector modules and the fiber optics interface attaching to the tomosynthesis system on the right.

Close-up view of the optical probes.

The TOBI system developed at MGH consists of two subsystems: the RF modulated imaging system and the CW system. The RF unit provides two laser wavelengths (685 and 830 nm) modulated at 70 MHz [10]. Optical switches (DiCon, Richmond, CA) multiplex the RF signal into 40 source locations and the light scattered through the tissue is collected by 8-RF avalanche photodiode detectors (APD, Hamamatsu, Bridgewater, NJ). The CW multiplexer unit (MUX) has three frequency-encoded lasers at 685, 810, and 830 nm. A fast Galvo scanner (Innovations in Optics Inc., Woburn, MA) is used to deliver the CW lasers to a maximum of 110 available locations on the optical probe with an average dwell time of 200 ms per location. A second frequency encoded CW system (TechEn Inc., Milford, MA) [22] provides 26 frequency encoded lasers, split equally between 685 and 830 nm, which are used for continuous monitoring of tissue changes at 26 source locations. The CW units share 32 APD (Hamamatsu C5460-01) detectors and the CW signals are subsequently demodulated for each channel and wavelength. For all MUX CW signals, de-multiplexing follows to separate the signals for each multiplexed source. The optical emission power exiting from the probe ranges from 10 to 20 mW for various RF and CW lasers. In addition to the optical system, a high-resolution linear encoder (Unimeasure Inc., Corvallis, OR) and four pressure transducers (Omigadyne Inc., OH) are mounted on the optical probes to give accurate readings of source-detector separation and applied pressure, respectively.

The TOBI optical probes (see Fig. 3) were specially designed to provide the capability of co-registration with 2-D mammography or tomosynthesis (more specifically, for GE units, though it could easily be adapted for other mammography machines). The anodized aluminum source probe covers a total area of 20 cm × 18 cm with fiber optics mounting locations over a 5-mm spacing grid. The metallic part of the probe can be securely attached to an optically transparent cassette modified from a standard mammography compression paddle. This cassette is fastened on the mammography machine. Similarly, the metallic detector probe can also be easily mounted on and detached from the moving compression arm using a push-button lock. The source/detector probes are inserted into the cassettes before taking the optical measurement and removed before taking X-ray scans. In both cases, the compression paddles retain their positions to ensure no movement of the target breast. Three groups of fiber optic bundles, mapped to different areas on the detector probe with partial overlaps, extend from the detector probe, terminated by metallic planar fiber connectors. These connectors can be interchangeably mounted to a female planar connector, aligned by mounting holes and magnet, and guide the transmitted light to RF and CW detector units. This design made it possible to conveniently switch optical scanning area based on mammographical views, such as cranio-caudal (CC), left mediolateral-oblique (LMLO), or right mediolateral-oblique (RMLO).

The tomosynthesis unit used in this study was developed by GE Healthcare, with the capability of collecting 15 projections within a 90° swing angle. The 3-D DBT images were reconstructed offline by a maximum-likelihood-type algorithm [23], [24]. The spatial resolution for a typical DBT image is 0.1 mm in the x and y axes and a z slice thickness of 1 mm.

B. Experiment Protocols

Two protocols were used in the clinical trial of this study. The operating steps for the protocols were as follows.

Protocol 1

1. Mount optical probes to X-ray compression paddles.
2. Compress patient breast to desired strength with compression paddles.
3. Perform optical data acquisition.
4. Remove optical probes from paddles while leaving patient's breast unmoved.
5. Perform DBT scan.
6. Release compression.
7. Repeat the above process for the contralateral breast (optional).

Protocol 2

1. Compress patient breast with compression paddles.
2. Perform DBT scan.
3. Mount optical probes to X-ray compression paddles while the patient breast remaining in compression.
4. Perform optical data acquisition.
5. Release compression.
6. Repeat the above process for the contralateral breast (optional).

In both cases, a solid calibration phantom measurement (optical only) is required before or after the patient data acquisition. For most of the experiments, more than one RF and CW source scan is completed within a 45 s duration. Typically the RF sources are scanned twice at 20 locations, while the MUX CW sources are scanned 7 times at 28 locations. The repeated measurements were averaged and used for the image reconstructions in this paper. The major difference between the two protocols is the temporal duration between the initial breast compression and the optical data acquisition. In Protocol 2, the elapsed time between the initial compression and optical data acquisition is roughly 1 or 2 min longer than in Protocol 1. Based on our study of breast compression induced tissue changes [25], we expect more hemodynamic change during the measurement period in Protocol 1 than in Protocol 2. However, tissue dynamics are not the major concern of this paper, a preliminary study on spatio-temporal image reconstruction of tissue dynamics can be found in Boverman et al. [26].

C. Data Analysis Procedure

Our data analysis procedures are outlined in the flow chart shown in Fig. 4. Most of these processing steps have been streamlined by software tools and require only minimum operator interference.

Data analysis flow chart.

The first step after retrieving the DBT image is to construct a coordinate mapping (or registration) between the DBT image voxels and the optical probe coordinates. The physical positions of curved breast boundaries are subsequently extracted from the registered DBT scan, from which the unstructured FEM mesh is generated with the simple algorithm described in Section II-D2. In most cases, not all of the optical sources/detectors were covered by the target breast. The stray light measurements from the uncovered sources/detectors therefore need to be removed from the raw data based on the registered DBT images and fiber locations.

The image reconstruction in this paper was done in two steps. In the first step, we estimated the bulk optical properties of the patient's breast as well as the mean source/detector coupling coefficients [21] (described in Section II-D). In the second step, we performed a full image reconstruction starting with the homogeneous initial guess from the output of the previous step. An indirect reconstruction scheme was used to obtain the hemoglobin concentrations and oxygen saturation (SO2) of the breast, meaning that we reconstructed the absorption and reduced scattering coefficients (μs′) of breast tissue at multiple wavelengths and then computed the hemoglobin concentration and SO2 based on HbR/HbO absorption spectra with assumed water and lipid concentrations [10], [27], [28].

D. Image Reconstruction Algorithm

In both the bulk optical property estimation and image reconstruction, we employed an iterative Gauss-Newton reconstruction approach. The forward problem is modelled by the finite element method (FEM) [14] under the diffusion approximation. The forward solutions at all valid sources and detectors are used to construct the Jacobian matrix using the adjoint method [29]-[31]. The absorption and reduced scattering coefficients of the target breast are iteratively updated by solving the normal equation. A major complication in this process is the simultaneous reconstruction of the source/detector coupling coefficients [21]. In this section, we outline the principal mathematical models used in the forward and inverse phases of the reconstruction paying especial attention to the coupling coefficient estimation.

1) Forward Model

A mathematical model for light propagations in human tissue is the radiative transport equation (RTE). In the regions where tissue is highly scattering, as is the case for breasts, RTE can be well approximated by a diffusion equation (DE). Assuming the light fluence rate is Φ, the diffusion equation for light propagation in the tissue can be written as [14]

−∇⋅D(r)∇Φ(r,t)+μa(r)Φ(r,t)+1c∂Φ(r,t)∂t=S0(r,t)

(1)

where μ_a_ is the absorption coefficient, D(r)=1∕3(μa(r)+μs′(r))≈1∕3μs′(r) is the diffusion coefficient and μ′s is the reduced scattering coefficient, c = _c_0/n is the speed of light in the medium, _c_0 is the speed of light in free-space, and n is the medium refraction index. _S_0(r,t) denotes the light source. For harmonic sources, with assumed time dependence exp(j_ω_t), the frequency-domain diffusion equation can be written as

−∇⋅D(r)∇Φ(r)+(μa(r)+jωc)Φ(r)=S0(r)

(2)

where ω is the angular frequency, Φ(r) is the phasor of Φ(r,t), and _S_0(r) is the phasor of the source S_0(r,t. Discretizing the continuous form by expanding Φ(r) on the 3-D forward mesh by Φ(r)=∑i=14Φiφi(r) along with expanding optical parameters on the reconstruction mesh by D(r)=∑k=14Dkϕk(r),μa(r)=∑k=14μakϕk(r), respectively, followed by integrating with weight function φ_j, we can get the FEM weak-form of (2) as

∑iΦi∑k(Dk〈ϕk,∇φi⋅∇φj〉−Dk〈ϕk∇φi,φj〉∂Ω+(μak+jωc)〈ϕkφi,φj〉)=〈S0,φj〉

(3)

where ϕ_k_ is the basis on the reconstruction mesh, φ_i_ and φ_j_ are the basis and weight functions on the forward mesh (assuming Galerkin's method [32]), respectively, the notation ⟨_f,g_⟩ denotes the intergration ∫Ω fgdr over volume Ω, and 〈f→,g〉∂Ω represents the surface integration ∫∂Ωgf→⋅dS→ On the boundaries ∂Ω.

On the boundaries, a partial-current boundary condition [33] was applied to correctly account for the refraction index discontinuity. This boundary condition can be expressed as [34]

where α is defined as

and _R_eff is the effective reflection coefficient [33]. Incorporating (4) into (3), we get the final form for the forward equation

∑iΦi(∑kDk〈ϕk,∇φi⋅∇φj〉+12α〈φi,φj〉∂Ω+∑k(μak+jωc)〈ϕkφi,φj〉)=〈S0,φj〉.

(6)

Assembling (6) for each element, we obtain a linear equation in form of Ax = b which can be solved by any linear equation solver.

2) Mesh Generation

A simple mesh generation approach was used to model the 3-D volume of the target breast. With this approach, we first create a uniform 3-D grid, and then split each cube based on a T5 rule (i.e., splitting a cube into five tetrahedrons) [35]. The resulting tetrahedral mesh is subsequently masked and truncated to fit the 3-D DBT images. Sometimes, a scaling along the _z_-axis is performed to match the probe separation in the experiment. To correctly consider the chest-wall area outside the DBT field-of-view, we also extend the DBT image toward the chest-wall direction for 2 cm before mesh generation. In this study, we used a dual-mesh scheme [36] in our reconstruction, i.e., a forward mesh for diffusion modeling and a separate reconstruction mesh to represent the optical properties. In particular, we create two sets of meshes based on the above process with different mesh densities, the higher one for the forward mesh, the lower one for the reconstruction mesh. The mesh files and the bilateral basis function mapping files are generated for each subject once and saved for all subsequent reconstructions.

3) Inverse Problem

The reconstruction of the tissue absorption and scattering properties is achieved by solving a nonlinear optimization problem, with the object function defined by

arg minμ_a_,D∥y−Φ(μ_a_,D)∥2 + λ∥(μ_a_,D)∥2

(7)

where y is the measurement vector and Φ(μa, D) is the data computed with the diffusion model (2) per given μa, and D; λ is a Tikhonov regularization parameter selected manually and ǁ·ǁ denotes the _L_2 norm or the sum of squares of log-amplitude/unwrapped phase differences [37]. The Gauss-Newton method iteratively updates the μa and D vectors by solving the following equation at each iteration:

(JkTJk+γkI+λI)(Δkμa,ΔkD)T=JkT(y−Φk(μa,D))

(8)

where △kμa and △k D are the property updates of absorption and diffusion coefficients, respectively, at the k_th iteration, γ_k is a Levenberg-Marquardt-type parameter determined by an empirical approach [38] at each iteration, I is an identity matrix. _J_k is the Jacobian matrix defined as

From (6), the submatrixes Jμa and JD can be computed by the adjoint method as [14], [29], [30], [39], [40]

Jμa((s,r),τ)=∂Φs,r∂μaτ=∑e∈Ωτ(MτeΦse)TΦ^re

(10)

JD((s,r),τ)=∂Φs,r∂Dτ=∑e∈Ωτ(HτeΦse)TΦ^re

(11)

where Φs,r represents the fluence rate measured at the _r_th detector under the illumination of the _s_th source, Σe∈Ωτ denotes the summation over all forward elements that fall in the immediate vicinity of the τth parameter node (in another word, where ϕτ is nonzero); the elements of matrixes Mτe and Hτe have forms of

mi,jτ=〈ϕτ(r)φi(r),φj(r)〉

(12)

hi,jτ=〈ϕτ(r),∇φi(r)⋅∇φj(r)〉

(13)

and vectors Φse and Φ^re are the Φ values at the four vertexes of the _e_th forward element under the illumination of the _s_th source and a unitary source at the _r_th detector (also referred to as an adjoint source), respectively.

In our work, using the log-amplitude and unwrapped-phase differences between y and Φ(μa, D) as the right-hand side of (8) was found to be superior in convergence [37] for bulk optical property estimation. The corresponding Jacobian matrix can be written as [40]

JνLP=1∥Φ∥2[(ℜ(Φ)ℜ(Jν)+ℑ(Φ)ℑ(Jν))+j(ℑ(Φ)ℜ(Jν)−ℜ(Φ)ℑ(Jν))]

(14)

where ν represents either μ_a_ or D; ℜ( ⋅ ) and ℑ( ⋅ ) are the real and imaginary parts, respectively, of a complex quantity. For CW data, we simply set the imaginary part of (14) to zero as the data only provide amplitude information.

Moreover, we applied a heuristic re-weighting scheme when solving (8). The weighted update equation can be written as

(JkTJk+γkΣ1+λΣ2)(Δkμa,ΔkD)T=JkT(y−Φk(μa,D))

(15)

where Σ1=diag(tr(JμaTJμa)1,tr(JDTJD)1) and Σ2 = diag({Σ_i_∣J i,j_∣2}j). The unitary vector 1 has the identical length of μa or D and Σ_i_∣_J i,_j_∣2 is the squared _L_2-norm of the _j_th column of the Jacobian matrix.

4) Data Calibration and Source/Detector Coupling Coefficient Estimation

Assuming a multiplicative model for the source-detector coupling coefficients [21], an appropriate calibration scheme is to divide the measurement of the target by that of a calibration phantom, and subsequently restore the calibration phantom data by multiplying the model prediction of the phantom. This can be written down as

ΦTcal=ΦTmeasΦPmeasΦPmodel

(16)

where ΦTmeas,ΦPmeas, and ΦPmodel denote the measured target data, measured, and model prediction of calibration phantom data, respectively. For diffuse optical tomography, the simple calibration scheme in (16) usually is not sufficient to remove all the systematic errors due to modeling approximations and experiment limitations. The uncalibrated systematic errors typically are manifested as artifacts aggregating around source/detector locations [21], [41]. In addition, TOBI does not use a liquid coupling bath as found in some DOT systems [42], [43]. As a result, the surface condition and contact between the calibration phantom and the breast measurement may vary significantly. This magnifies the source/detector artifacts. Thus, it is critically important to choose an appropriate model for the coupling coefficients to obtain meaningful images in our study.

As proposed by Boas et al. [21] and further explored by Oh et al. [44], Tarvainen et al. [45], and Schweiger et al. [46], we assumed a multiplicative measurement model with unknown source/detector coupling coefficients

Φi,j=si×dj×Φi,jmodel(μa,D)

(17)

where si and dj are the multiplicative coupling coefficient for the _i_th source and _j_th detector, respectively. The image reconstruction problem involves the simultaneous determination of coupling coefficients, s = {si} and d = {d j} , along with the estimation of optical properties, μa and D. The normal equation is now reformulated as

(JkTJk+γkI+λI)(Δkμa,ΔkD,Δks,Δkd)T=JkT(y−diag(s⊗d)Φk(μa,D))

(18)

where

and ⊗ denotes the Kronecker product. From (17), it is not difficult to show

Simultaneously estimating source/detector coupling coefficients and optical images greatly reduces the previously mentioned artifacts [21]. However, a side-effect of this scheme is the reduction of target contrast, particularly when measurements were collected within limited projection angles. In these cases, a factor of 2 or higher of contrast reduction can be found in simulations and phantom data reconstruction. This is not surprising since (18) is effectively estimating three images simultaneously: the 3-D optical property distributions, the source coupling coefficient distribution at the source plane and analogously for detectors. Without properly constraining the unknowns, the contrast from one image could be easily “stolen” by other images. In order to get images with appropriate optical contrast, a priori knowledge about s and d should be applied when solving (18). If the surface condition and fiber contact are fairly close between calibration phantom and patient, we can reasonably assume that SD coefficients are likely to have norm 1 on average. Mathematically this constraint can be represented by adding an additional Tikhonov regularization term into the object function as

argminμa,D,s,d∥y−diag(s⊗d)Φ(μa,D)∥2+λ∥(μa,D)∥2+α∥(s,d)T−1∥2

(22)

where α is an additional regularization parameter controlling the strength of this constraint. Meanwhile, we noticed that applying a spatial high-pass filter to the SD distribution is useful to reduce the crosstalk between the optical images and SD coefficients, for the fact that our optical images contain mostly the low spatial-frequency components. This leads to a third regularization term in the object function

argminμa,D,s,d∥y−diag(s⊗d)Φ(μa,D)∥2+λ∥(μa,D)∥2+α∥(s,d)T−1∥2+β∥Γ(s,d)∥2

(23)

where matrix Γ is a discrete integration operator. For the cases where the number of measurements is less than the number of unknowns, i.e., the underdetermined cases, the presence of the third term in (23) could lead to a increase in computational time. A low-cost alternative is to perform a high-pass filter on SD distributions as a post process in each iteration. The filtered source coefficient can be expressed as

where the index j loops over the neighbors of the _i_th source from a tesselated source location mesh. A similar expression can be used for the detector coefficients. Although the filtering scheme in (24) is not as rigorous as the regularization approach in (23), numerical simulations show that the results from the filtering approach are reasonably close to the true values. Equation (24) was used for all reconstructions in Section III.

5) Bulk Property Estimation and Image Reconstruction

With the derived reconstruction scheme, we estimated the breast bulk optical properties from the RF measurements. This is done by regarding the whole reconstruction mesh as a single set of (μa,D)pairs. In this case, the Jacobian matrix is equal to the sum of columns for all parameter nodes. The SD unknowns are also updated simultaneously in the bulk property estimation to absorb the systematic errors.

The estimated bulk optical properties were used as homogeneous initial guess for the image reconstruction stage. However, the SD estimated from the bulk estimation can not be used for the second stage. This is because the SD coefficients in the bulk estimation have maximum coupling to the optical images since they are the only degree-of-freedom used to represent spatial variability. However, averages of source and detector coefficients can be used to initialize the SD coefficients in the image reconstruction.

III. RESULTS

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], [47]-[49].

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 Fig. 5. 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 Fig. 5(a) and (b). 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 Fig. 6(a). 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 Fig. 6 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 [Fig. 6(b)].

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 and μs′=5.48cm−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 and μs′=6.00cm−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],[51]. 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 Fig. 7. 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.

Reconstructed image slices for a balloon phantom filled with intralipid solution: (a)μ_a_ (cm-1) image and (b) μs′(cm−1) image. Both images were extracted at the horizontal plane 3 cm from the bottom of the phantom. Maximum perturbations for the absorption is ±2.1% and that for the scattering image is ±1.9%.

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 Fig. 6(c). 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 Fig. 8. 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.

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 overlayed on the images.

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 (μs′) at 830 nm are listed in Table I. 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], [47]-[49]. However, due to the applied mammographic compression, removal of blood is anticipated from our previous studies [25] and the values in Table I 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 μs′

Twelve DBT image slices (six for LMLO view, six for RMLO view) are shown in Fig. 9 (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 Fig. 9(b)-(f) and Fig. 9(i)-(j)].

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 dense cases. All images are slices extracted along a horizontal plane from the reconstructed 3-D HbT maps. The thick dashed line denotes the edge of optical source/detector coverage. The region between the dashed line and the nipple has virtually no sensitivity to the measurements. The white dashed line on the DBT slice and the black solid thick line on the HbT images denote the boundaries of the chest-wall muscle regions which do not show up for all the patients. In all figures, the color scales were set between 80%-120% of the respective bulk HbT value, except for (g) where 70%-130% was used.

Reconstructed HbT and DBT images for 12 representative healthy breasts are shown in Fig. 9. 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], [55], 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 Fig. 9. For breasts with very little glandular content, i.e., Fig. 9(a)-(c), 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 Fig. 9(d)-(f), 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 Fig. 9(i), 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 Fig. 9(i). The HbT decrease in these regions ranges from about 5% to about 20%.

To generalize the findings from Fig. 9, 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.

The reconstructed μs′ distributions are similar to those of HbT for each patient. Representative images of the μs′ from two breasts, i.e., Fig. 9(i) and (d), are shown in Fig. 10(a) and (c), respectively, where the glandular regions demonstrate higher scattering values. The corresponding SO2 images of the two subjects [Fig. 10(b) and (d)] seem to present lower values within the fibroglandular segmentations, which suggest higher metabolic activities in the fibroglandular tissue.

Image slices for (a) μs′(cm−1) at 830 nm, (b) oxygen-saturation (SO2) corresponding to case (i) in Fig. 10, and (c) μs′(cm−1) at 830 nm, (d) SO2 corresponding to case (d) in Fig. 10. The markings on the images are defined similarly as in Fig. 9.

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 Fig. 11(a). 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 and y 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 Fig. 11(b) and (c). 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], [58]. 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.

IV. CONCLUSION

We have demonstrated a working prototype of a combined X-ray and optical breast imaging system and its successful application to imaging healthy breasts in a clinical setting. The capability of performing optical and X-ray scans in a spatially co-registered manner allows us to interpret the optical images with greater confidence. The reconstructed HbT images of breasts were reasonably consistent with the corresponding tomosynthesis images. The contours of the chest-wall muscle, fibroglandular and fatty tissue regions in the optical images generally agree with the X-ray structural information. An ROI analysis has demonstrated significant differences in HbT between these tissue types. Particularly, the average contrast between the fibroglandular and adipose tissue is 1.4:1. For breasts with very small and very large glandular regions, we observed a slight decrease in the central part of the HbT images. We hypothesized that this is due to the nonuniformity of the applied mammography compression and the induced hemoglobin redistribution. A mechanical simulation of breast compression qualitatively strengthened this hypothesis.

Despite the fact that the optical images have low spatial resolution, interpreting the features in the images is fairly straightforward with the guidance of co-registered X-ray scans. As one can see, the combined optical/X-ray breast imaging system presents additional diagnostic information in an incremental and co-registered manner. This could possibly lead to a lower barrier for acceptance by radiologists.

The blood-pressure coupling effect revealed in this paper and some previous work [25] opens a new horizon for breast cancer detection. On the one hand, the precisely controlled external compression provide a simple exogenous contrast mechanism for hemoglobin concentration, where both the tissue mechanical and vascular characteristics were encoded and can be detected by optical measurement. It may be possible to fine-tune the system so that tumor abnormalities in both the mechanical map [59] and the hemoglobin map enhance each other, resulting in a more efficient detection. On the other hand, the tissue metabolism and the visco-elasticity [56] add a new dimension, spatial-temporal dynamics, to cancer identification where dynamic compression and real-time optical measurements are simultaneously available. To study the temporal signature of the tissue under compression is also a possible path to improve the imaging specificity.

Our next goal is to study the system's performance for breasts with lesions, as this is obviously critical for proving the system's clinical significance. As we concluded from the healthy subject study, mammographic compression affects the image contrast in some cases. It is natural to ask whether the compression will affect tumor detection. To answer this question, we will reconstruct images for tumor-bearing breasts. Using the ROIs derived from X-ray images, we can calculate the optical contrast and the statistical significance between the tumor and healthy tissues, and then compare with the published optical contrast from DOT systems without compression [61]. Tumor segmentation information will be used in the reconstruction as a spatial prior [62] to reduce the underestimation for small targets as shown in our phantom reconstruction. In the meantime, our breast compression study revealed interesting tissue dynamics such as pressure relaxation, oxygenation and flow changes associated with breast compression [25]. Our TOBI measurements using various protocols (Section II-B) also indicated an up to 5% monotonic signal change during the experiment (the variations in Protocol 1 are slightly greater than those in Protocol 2). Although these variations are not large enough to alter the results of the static image reconstruction, incorporating dynamic models into the analysis [26] will give us opportunities to study the temporal signature of the tissue dynamics and provide new dimension for differentiating cancers from healthy tissue in the early stage.

Acknowledgments

This work was supported in part by National Institutes of Health (NIH) through the National Cancer Institute under Grant R01-CA97305-01 and Grant U54-CA105480 (NTROI) and in part by the Center for Subsurface Sensing and Imaging System (CenSSIS) under the Engineering Research Centers Program of the National Science Foundation Award EEC-9986821.

REFERENCES

[2] Tabar L, Yen M-F, Vitak H-HTCB, Smith RA, Duffy SW. Mammography service screening and mortality in breast cancer patients: 20-year follow-up before and after introduction of screening. Lancet. 2003;361:1405–1410. [PubMed] [Google Scholar]

[3] Hebden JC, Arridge SR, Delpy DT. Optical imaging in medicine: I. Experimental techniques. Phys. Med. Biol. 1997;42:825–840. [PubMed] [Google Scholar]

[4] Boas DA, Brooks DH, Miller EL, DiMarzio CA, Kilmer M, Gaudette RJ, Zhang Q. Imaging the body with diffuse optical tomography. IEEE Signal Process. Mag. 2001 Nov.18(6):57–75. [Google Scholar]

[5] Gibson AP, Hebden JC, Arridge SR. Recent advances in diffuse optical imaging. Phys. Med. Biol. 2005;50:R1–R43. [PubMed] [Google Scholar]

[6] Zhu Q, Chen NG, Kurtzman SH. Imaging tumor angiogenesis using combined near infrared diffusive light and ultrasound. Opt. Lett. 2003;28:337–339. [PMC free article] [PubMed] [Google Scholar]

[7] Pogue BW, Poplack SP, McBride TO, Wells WA, Osterman KS, Osterberg UL, Paulsen KD. Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: Pilot results in the breast. Radiology. 2001;218:261–266. [PubMed] [Google Scholar]

[8] Shah N, Gibbs J, Wolverton D, Cerussi A, Hylton N, Tromberg BJ. Combined diffuse optical spectroscopy and contrast-enhanced MRI for monitoring breast cancer neoadjuvant chemotherapy: A case study. J. Biomed. Opt. 2005;10(5):051503–051503. [PubMed] [Google Scholar]

[9] Choe R, Corlu A, Lee K, Durduran T, Konecky SD, Grosicka-Koptyra M, Arridge SR, Czerniecki BJ, Fraker DL, DeMichele A, Chance B, Rosen MA, Yodh AG. Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: A case study with comparison to MRI. Med. Phys. 2005;32:1128–1139. [PubMed] [Google Scholar]

[10] Zhang Q, Brukilacchio TJ, Li A, Stott JJ, Chaves T, Hillman E, Wu T, Chorlton M, Rafferty E, Moore RH, Kopans DB, Boas DA. Coregistered tomographic X-ray and optical breast imaging: Initial results. J. Biomed. Opt. 2005;10:024033–024033. [PubMed] [Google Scholar]

[11] Brooksby B, Jiang S, Dehghani H, Pogue BW, Paulsen KD, Weaver J, Kogel C, Poplack SP. Combining near infrared tomography and magnetic resonance imaging to study in vivo breast tissue: Implementation of a laplacian-type regularization to incorporate magentic reasonance structure. J. Biomed. Opt. 2005;10:0515041–10. [PubMed] [Google Scholar]

[12] Cerussi A, Shah N, Hsiang D, Durkin A, Butler J, Tromberg BJ. In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy. J. Biomed. Opt. 2006;11 [PubMed] [Google Scholar]

[13] Klifa CS, Hattangadi J, Li A, Tromberg BJ, Hylton NM, Park C. Combination of magnetic resonance imaging and diffuse optical spectroscopy to predict radiation response in the breast: An exploratory pilot study. Breast Cancer Res. Treatment. 2006;100:S202–S202. [Google Scholar]

[14] Arridge SR. Optical tomography in medical imaging. Inv. Problems. 1999;5:R41–R93. [Google Scholar]

[15] Teicher BA. Physiologic mechanisms of therapeutic resistance: Blood flow and hypoxia. Hematol. Oncol. Clin. North Am. 1995;9:475–506. [PubMed] [Google Scholar]

[16] Vaupel P, Kallinowski F, Okunieff P. Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: A review. Cancer Res. 1989;49:6449–6465. [PubMed] [Google Scholar]

[17] Czernin J, Schelbert HR. PET/CT in cancer patient management. J. Nucl. Med. 2007;48 [PubMed] [Google Scholar]

[18] Fang Q, Boas DA, Boverman G, Zhang Q, Kauffman T. Nonlinear image reconstruction algorithm for diffuse optical tomography using iterative block solver and automatic mesh generation from tomosynthesis images. Proce. SPIE. 2006 Feb.6081 [Google Scholar]

[19] Boyse WE, Seidl AA. A block QMR method for computing multiple simultaneous solutions to complex symmetric systems. SIAM J. Sci. Comput. 1996;17:263–274. [Google Scholar]

[20] Joachimowicz N, Pichot C, Hugonin JP. Inverse scattering: An iterative numerical method for electromagnetic imaging. IEEE Trans. Antennas Propag. 1991 Dec.39(12):1742–1752. [Google Scholar]

[21] Boas D, Gaudette T, Arridge S. Simultaneous imaging and optode calibration with diffuse optical tomography. Opt. Express. 2001;8:263–270. [PubMed] [Google Scholar]

[22] Joseph DK, Huppert TJ, Franceschini MA, Boas DA. Diffuse optical tomography system to image brain activation with improved spatial resolution and validation with functional magnetic resonance imaging. Appl. Opt. 2006;45:8142–8151. [PubMed] [Google Scholar]

[23] Wu T, Stewart A, Stanton M, McCauley T, Phillips W, Kopans DB, Moore RH, Eberhard JW, Opsahl-Ong B, Niklason L, Williams MB. Tomographic mammography using a limited number of low-dose conebeam projection images. Med. Phys. 2003;30:365–380. [PubMed] [Google Scholar]

[24] Wu T, Zhang J, Moore R, Rafferty E, Kopans D, Meleis W, Kaeli D. Digital tomosynthesis mammography using a parallel maximum-likelihood reconstruction method. Proc. SPIE. 2004;5368:1–11. [Google Scholar]

[25] Carp SA, Kauffman T, Fang Q, Rafferty E, Moore R, Kopans D, Boas D. Compression-induced changes in the physiological state of the breast as observed through frequency domain photon migration measurements. J. Biomed. Opt. 2006;11(6):064016–064016. [PubMed] [Google Scholar]

[26] Boverman G, Fang Q, Carp SA, Miller EL, Brooks DH, Selb J, Moore RH, Kopans DB, Boas DA. Spatio-temporal imaging of the hemoglobin in the compressed breast with diffuse optical tomography. Phys. Med. Biol. 2007;52:3619–3641. [PubMed] [Google Scholar]

[27] Paulsen KD, Jiang H. Spatially varying optical property reconstruction using a finite element diffusion equation approximation. Med. Phys. 1995;22:691–701. [PubMed] [Google Scholar]

[28] Arridge SR, Schweiger M. Image reconstruction in optical tomography. Philos. Trans. R. Soc. London Ser. B. 1997;352:717–726. [PMC free article] [PubMed] [Google Scholar]

[29] Vogel CR. Computational Methods for Inverse Problems. SIAM; Philadelphia, PA: 2002. [Google Scholar]

[30] Liauh C-T, Hills RG, Roemer RB. Comparison of the adjoint and influence coefficient methods for solving the inverse hyoperthermia problem. J. Biomechan. Eng. 1993;115:63–71. [PubMed] [Google Scholar]

[31] Fang Q, Meaney PM, Geimer SD, Paulsen KD, Streltsov AV. Microwave image reconstruction from 3-D fields coupled to 2-D parameter estimation. IEEE Trans. Med. Imag. 2004 Apr.23(4):475–484. [PubMed] [Google Scholar]

[32] Zienkiewicz OC, Taylor RL. Finite Element Method-Basic Formulation and Linear Problems. vol. 1. McGraw-Hill; New York: 1989. [Google Scholar]

[33] Haskell RC, Svaasand LO, Tsay T, Feng T, McAdams MS, Tromberg BJ. Boundary conditions for the diffusion equation in radiative transfer. J. Opt. Soc. Am. A. 1994;11:2727–2741. [PubMed] [Google Scholar]

[34] Arridge SR, Dehghani H, Schweiger M, Okada E. The finite element model for the propagation of light in scattering media: A direct method for domains with nonscattering regions. Med. Phys. 2000 Jan.27:252–264. [PubMed] [Google Scholar]

[35] Fleischmann P, Haindl B, Kosik R, Selberherr S. Investigation of a mesh criterion for three-dimensional finite element diffusion simulation. Proc. Int. Conf. Simulation Semiconductor Processes Devices (SISPAD '99) 1999:71–74. [Google Scholar]

[36] Paulsen KD, Meaney PM, Moskowitz MJ, Sullivan JM., Jr. A dual mesh scheme for finite element based reconstruction algorithms. IEEE Trans. Med. Imag. 1995 Sep.14(3):504–514. [PubMed] [Google Scholar]

[37] Meaney PM, Paulsen KD, Pogue BW, Miga MI. Microwave image reconstruction utilizing log-magnitude and unwrapped phase to improve high-contrast object recovery. IEEE Trans. Med. Imag. 2001 Feb.20(2):104–116. [PubMed] [Google Scholar]

[38] Hugonin JP, Joachimowicz N, Pichot C. Inverse Methods in Action. Springer-Verlag; Berlin, Germany: 1990. [Google Scholar]

[39] Chavent G, Lemmonier P. Identification de la non linéaritéd'uneéquation parabolique quasilinéaire. Appl. Math. Optimizat. 1974;1:121–162. [Google Scholar]

[40] Fang Q. Ph.D. dissertation. Dartmouth College; Hanover, NH: 2004. Computational methods for microwave medical imaging. [Google Scholar]

[41] Li C, Jiang H. A calibration method in diffuse optical tomography. J. Opt. A: Pure Appl. Opt. 2004;6:844–852. [Google Scholar]

[42] Durduran T, Choe R, Culver JP, Zubkov L, Holboke MJ, Giammarco J, Chance B, Yodh AG. Bulk optical properties of healthy female breast tissue. Phys. Med. Biol. 2002;47:2847–2861. [PubMed] [Google Scholar]

[43] Nielsen T, Brendel B, Koehler T, Ziegler R, Ziegler A, Bakker L, Beek M. v., Mark M. v. d., Voort M. v. d., Harbers R, Licha K, Pessel M, Schippers F, Meeuwse JP, Feuerabend A, Pijkeren D. v., Deckers S. Image reconstruction and evaluation of system performance for optical fluorescence tomography. Proc. SPIE. 2007;6431:643108.1–643108.10. [Google Scholar]

[44] Oh S, Milstein AB, Millane RP, Bouman CA, Webb KJ. Source detector calibration in three-dimensional bayesian optical diffusion tomography. J. Opt. Soc. Am. A. 2002;19:1983–1993. [PubMed] [Google Scholar]

[45] Tarvainen T, Kolehmainen V, Vauhkonen M, Vanne A, Gibson AP, Schweiger M, Arridge SR, Kaipio JP. Computational calibration method for optical tomography. Appl. Opt. 2005;44:1879–1888. [PubMed] [Google Scholar]

[46] Schweiger M, Nissilä I, Boas DA, Arridge SR. Image reconstruction in optical tomography in the presence of coupling errors. Appl. Opt. 2007;46:2743–2756. [PubMed] [Google Scholar]

[47] Cerussi AE, Berger AJ, Bevilacqua F, Shah N, Jakubowski D, Butler J, Holcombe RF, Tromberg BJ. Sources of absorption and scattering contrast for near-infrared optical mammography. Acad. Radiol. 2001;8:209–10. [PubMed] [Google Scholar]

[48] Srinivasan S, Pogue BW, Jiang S, Dehghani H, Kogel C, Soho S, Gibson JJ, Tosteson TD, Poplack SP, Paulsen KD. In vivo hemoglobin and water concentration, oxygen saturation, and scattering estimates from near-infrared breast tomography using spectral reconstruction. Acad. Radiol. 2006;13:195–202. [PubMed] [Google Scholar]

[49] Brooksby B, Pogue BW, Jiang S, Dehghani H, Srinivasan S, Kogel C, Tosteson TD, Weaver J, Poplack SP, Paulsen KD. Imaging breast adipose and fibroglandular tissue molecular signatures using hybrid MRI-guided near-infrared spectral tomography. Proc. Natl. Acad. Sci. USA. 2006;103:8828–8833. [PMC free article] [PubMed] [Google Scholar]

[50] Pifferi A, Torricelli A, Bassi A, Taroni P, Cubeddu R, Wabnitz H, Grosenick D, Moller M, Macdonald R, Swartling J. Performance assessment of photon migration instruments: The medphot protocol. Appl. Opt. 2005;44:2104–2114. [PubMed] [Google Scholar]

[51] Selb J, Joseph DK, Boas DA. Time-gated optical system for depth-resolved functional brain imaging. J. Biomed. Opt. 2006;11:044008–044008. [PubMed] [Google Scholar]

[52] Yalavarthy PK, Pogue BW, Dehghani H, Paulsen KD. Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography. Med. Phys. 2007 Jun.34:2085–2098. [PubMed] [Google Scholar]

[53] Yoo TS, Ackerman MJ, Lorensen WE, Schroeder W, Cha-lana V, Aylward S, Metaxes D, Whitaker R. Engineering and algorithm design for an image processing API: A technical report on ITK—The insight toolkit. Proc. Medicine Meets Virtual Reality. 2002:586–592. [PubMed] [Google Scholar]

[54] American College of Radiology (ACR) Breast Imaging Reporting and Data System Atlas (BI-RADS Atlas) Amer. College Radiol.; Reston, VA: 2003. [Google Scholar]

[55] Shah N, Cerussi A, Eker C, Espinoza J, Butler J, Fishkin J, Hornung R, Tromberg B. Noninvasive functional optical spectroscopy of human breast tissue. PNAS. 2001;98:4420–4425. [PMC free article] [PubMed] [Google Scholar]

[56] Fung YC. Biomechanics: Mechanical Properties of Living Tissues. 2nd ed. Springer Verlag; New York: 1993. [Google Scholar]

[57] Berry GP, Bamber JC, Miller NR, Barbone PE, Bush NL, Armstrong CG. Towards an acoustic model-based poroelastic imaging method: II. Experimental investigation. Ultrasound Med. Biol. 2006;32:1869–1885. [PubMed] [Google Scholar]

[58] Kerdok AE, Ottensmeyerc MP, Howe RD. Effects of per-fusion on the viscoelastic characteristics of liver. J. Biomechan. 2006;39:2221–2231. [PubMed] [Google Scholar]

[59] Wellman PS, Howe RD. The mechanical properties of breast tissues in compression. Harvard BioRobotics Lab. Tech. Rep. 1999;99003 [Google Scholar]

[60] Jain RK. Transport of molecules in the tumor interstitium: A review. Cancer Res. 1987 Jun.47:3039–3051. [PubMed] [Google Scholar]

[61] Poplack SP, Tosteson TD, Wells WA, Pogue BW, Meaney PM, Hartov A, Kogel CA, Soho SK, Gibson JJ, Paulsen KD. Electromagnetic breast imaging: Results of a pilot study in women with abnormal mammograms. Radiology. 2007 May;243:350–359. [PubMed] [Google Scholar]

[62] Li A, Miller EL, Kilmer ME, Brukilacchio TJ, Chaves T, Stott J, Zhang Q, Wu T, Chorlton M, Moore RH, Kopans DB, Boas DA. Tomographic optical breast imaging guided by three-dimensional mammography. Appl. Opt. 2003;42:5181–5190. [PubMed] [Google Scholar]