Open Access
Add to BibSonomyAbstract
Multi-spectral near-infrared diffuse optical tomography (DOT) is capable of providing functional tissue assessment that can complement structural mammographic images for more comprehensive breast cancer diagnosis. To take full advantage of the readily available sub-millimeter resolution structural information in a multi-modal imaging setting, an efficient x-ray/optical joint image reconstruction model has been proposed previously to utilize anatomical information from a mammogram as a structural prior. In this work, we develop a complex digital breast phantom (available at http://openjd.sf.net/digibreast) based on direct measurements of fibroglandular tissue volume fractions using dual-energy mammographic imaging of a human breast. We also extend our prior-guided reconstruction algorithm to facilitate the recovery of breast tumors, and perform a series of simulation-based studies to systematically evaluate the impact of lesion sizes and contrasts, tissue background, mesh resolution, inaccurate priors, and regularization parameters, on the recovery of breast tumors using multi-modal DOT/x-ray measurements. Our studies reveal that the optical property estimation error can be reduced by half by utilizing structural priors; the minimum detectable tumor size can also be reduced by half when prior knowledge regarding the tumor location is provided. Moreover, our algorithm is shown to be robust to false priors on tumor location.
© 2015 Optical Society of America
1. Introduction
Diffuse optical tomography (DOT) is a functional imaging modality used for recovery of deep tissue physiological parameters, such as tissue hemoglobin concentrations and oxygenation, by using only non-ionizing near-infrared (NIR) radiation [1]. Due to the low overall tissue absorption within the “optical window” spanning from around 600 to1000 nm, NIR light is capable of penetrating through over a dozen centimeters of human tissue [2]. In particular, DOT has been extensively studied for breast imaging applications [3,4], including imaging breast cancer during diagnosis [5–10] and monitoring treatment response during neoadjuvant chemotherapy [11,12]. This is primarily due to the fact that the breast has good accessibility for making optical measurements and has overall low optical absorption due to the high percentage of fatty tissue. However, because DOT image reconstructions suffer from ill-posedness, i.e. intrinsic sensitivity to measurement noise and model errors [13,14], DOT reconstructions generally exhibit low spatial resolution [15]. For the tissue thicknesses that are commonly found in breast imaging, the minimum feature size that can be reliably resolved by DOT has been reported to be typically between 7 to 10 mm [16,17]. The low spatial resolution of DOT greatly limits its more widespread adoption in the clinic.
Thus investigating new methods to improve DOT spatial resolution could have significant implications for increasing the utility of this technique for clinical applications where image spatial details are critical. Increasing the density of measurements may be one way to help to improve the resulting resolution. Examples of this approach include the recent development of high-density optical probes [16] and camera-based detection approaches [17–19]. However, the ill-posed nature of DOT image reconstruction limits the effectiveness of using more and more measurements – the more data that are included, the less effective the increase is for improving image resolution. One alternative approach is to use structural priors [20]. By incorporating spatial information from another modality, preferably a structure-oriented image [21–23], the optical reconstructions can be “guided” by the tissue anatomy and may result in more accurate spatial details. While such structural guidance may be achieved during image interpretation, where the both image sets are simply shown to the clinician simultaneously and interpreted together, using the two imaging data sets jointly during image reconstruction represents a step forward in full exploitation of the imaging information. In DOT breast imaging, the latter idea has been explored by using binary-segmented structural-regions to either reduce the number of unknowns in the image space (known as the “hard” prior [24]), or to impose a region-based smoothing (known as the “soft” prior [25–27]).
Together with several other groups [28,29], our group has further advanced the prior-guided reconstruction technique by proposing a generalized structural-prior guided reconstruction [30]. Instead of forcing the background tissue structure into fixed, piecewise-constant regions, or even pre-determined linear combinations of such regions [29], we apply a “fuzzy” segmentation of the tissue anatomy, represented by a probability (or volume fraction) map for each tissue composition. This approach preserves the fine spatial details from the structural images, and has shown robustness in processing clinical measurements [30]. Nonetheless, a rigorous validation of this approach has not been reported, partly due to the lack of clinically realistic optical breast phantoms. Although our in vivo study from a clinical population has anecdotally demonstrated image resolution enhancement using this approach [30], without knowing the ground-truth of the optical properties in a real breast, it is quite difficult to systematically quantify the performance of the algorithm, especially when we also consider the inaccuracy of the structural priors themselves.
In recent years, the development of photon-counting spectral imaging in the field of mammography [31] has enabled the direct quantification of two-dimensional (2D) breast tissue compositions, i.e. adipose and fibroglandular tissues, from one standard mammography acquisition [32]. Such compositional information is precisely what is needed in our prior-guided optical reconstruction algorithm. A clinically measured tissue compositional map can also serve as the “ground-truth” to validate our algorithm using simulations. While a 2D tissue compositional map remains different from a real breast anatomy, it offers significantly more spatial details of the breast tissue structures compared to the uniform or binary-segmented breast models used in the literature [20–29], and therefore represents a significant step forward towards validating parallel-plate DOT breast imaging [7] in more realistic conditions.
The goal of this paper is to systematically characterize the performance and optimize the key parameters of our previously proposed prior-guided DOT reconstruction algorithm, in particular for better recovery of small and low-contrast tumors. To do this, we first construct a digital breast phantom using clinically derived 3D tissue compositional maps, and then run a series of simulations to test the impact of 1) structural image fuzzy segmentation algorithms, 2) regularization, 3) mesh density, 4) tumor size, and 5) tumor-to-background contrast. We hope that the findings from this work will be generalizable to other multi-modal DOT imaging systems such as those combining DOT with MRI [21,22] and ultrasound [8].
In the following sections, we first describe the key methods of this analysis, including the creation of the digital breast phantom, forward modeling and measurement simulation, along with a number of extensions to our image reconstruction algorithm; then we outline the various benchmarks performed in this simulation-based study. The outcomes from these benchmarks are summarized in the Results section. Finally, our key findings regarding the algorithm performance and parameter optimization are summarized in the Discussions section.
2. Methods
2.1 Digital breast phantom derived from clinical dual-energy x-ray scans
A two-dimensional digital breast mammogram was acquired using a Philips MicroDose SI mammography system, a scanning multi-slit system with photon-counting silicon strip detectors. A threshold in the x-ray detector electronics sorts detected photons into two bins according to photon energy [31]. By calibrating against a known pair of materials, it is possible to map the photon counts in both detection bins to any other known pair of materials [33]. As a result, in addition to the 2D digital mammogram, this MicroDose SI system is capable of measuring the breast tissue compositions [32] by calibrating and validating the counts in the two bins on two known breast density phantom materials and then mapping to published attenuation of adipose and fibroglandular tissue [34].
Using this technique, we have obtained a 2D fibroglandular tissue volume fraction map, referred to as Cf, from a normal breast. The Cf image, shown in Fig. 1(b) , has the same size and view as the 2D mammogram (Fig. 1(a)), with pixel values between 0 and 1 that represent the percentage volume fractions of the fibroglandular tissue. Along with Cf, the scanner also provides the thickness of the mammographically compressed breast, denoted here by T, as a function of pixel location, making it possible to compute the total breast volume. The 2D cranial-caudal (CC) view mammogram and the corresponding Cf image are shown in Figs. 1(a) and 1(b) respectively. Note that the skin region that surrounds the outer contour is removed from the displayed image to minimize the bias due to skin attenuation. The compression thickness of the breast is 30 mm; the total area in this view is 161.8 cm2, and the total volumetric percentage of the fibroglandular tissue is 24%, calculated based on both Cf and T.
Fig. 1 A 3D digital breast phantom derived from a dual-energy x-ray mammography scan: (a) 2D mammogram, (b) fibroglandular volume fraction, (c) 3D uniform breast mesh (color showing total hemoglobin concentrations, inset showing the vertical profile), and (d,e) cross-sectional views of locally-refined meshes along with optical fiber locations. Two simulated tumor locations, Ωa and Ωf, are indicated with arrows.
A 3D digital breast phantom was subsequently constructed using the Cf image and the thickness information. We first down-sampled Cf to a pixel size of 1 by 1 mm, and then expanded the down-sampled Cf into a 3D cylindrical structure by stacking repeated slices along the z-axis with an assumed slice thickness of 1 mm. This quasi-3D Cf profile was subsequently masked by a 3D breast shape profile derived by mirroring half of the breast thickness map, T/2, across the middle slice of the stack (shown as the inset in Fig. 1(c)). To minimize the optical modeling errors caused by volume truncation, we also add a 2 cm wide slab to the volume towards the chest-wall by replicating the Cf boundary pixels (the extended region was excluded from the error calculations due to the absence of optical fiber coverage).
Based on this 3D fibroglandular volume fraction profile, a pair of tetrahedral meshes, a finer one for solving the optical forward problem and a coarser one for solving the inversion, were generated using a MATLAB-based meshing toolbox “iso2mesh” [35]. Briefly, for either mesh, we first extracted a triangular surface from the volume and then smoothed this surface by applying a low-pass filter [36]. We then populated the tetrahedral elements inside the space bounded by the smoothed surface. The resulting forward mesh is shown in Fig. 1(c).
To quantify the performance of our algorithm, we test two different tumor-to-background contrasts in this study: 1) tumor is located inside the adipose tissue (Ωa in Fig. 1(d)), and 2) tumor is located inside the fibroglandular tissue (Ωf in Fig. 1(e)). In addition, to study the impact of mesh density on the minimum tumor size that can be recovered, we also construct a set of refined meshes using iso2mesh by increasing the mesh density within a 2 cm diameter sphere centered at either candidate tumor location (Ωa or Ωf). A summary of the forward and reconstruction meshes with the uniform and refined (tumor region Ωa) density settings is shown in Table 1 . Cross-sectional views of the refined mesh are shown in Figs. 1(d) and 1(e).
Table 1. Statistics of the forward and reconstruction meshes generated from the digital breast phantom
Using the x-ray derived 3D fibroglandular volume fraction map and the breast meshes, we then defined the optical properties at each node in the forward mesh by applying the following formula (thus implicitly assuming that the primary breast tissue constituents are fibroglandular and adipose tissue only):
where µ fib = {HbO, HbR, μs’(690 nm), μs’(830 nm)}fib denotes the set of optical properties, including the oxy-hemoglobin concentration (HbO), deoxy-hemoglobin concentration (HbR), and the reduced scattering coefficients (μs’) at 690 and 830 nm, for the fibroglandular tissue, and µ adi the same for the adipose tissue. Ca (r) is the adipose tissue volume fraction at location r, computed as 1- Cf (r). The optical property values for µ adi and µ fib were determined from averages of the recovered adipose and fibroglandular optical properties, respectively, using the clinical breast measurements (N = 189) reported in [9], and are reported in Table 2 .Table 2. Optical properties for various breast tissues
To test the algorithm performance for tumor recovery, we artificially added a 3D breast tumor to the above normal breast phantom. The tumor was characterized by a tumor tissue volume fraction map in the shape of a 3D Gaussian-sphere, represented by
where ||∙|| is the L-2 norm; r 0 is the centroid of the tumor; σ is the standard deviation of the Gaussian sphere. The effective diameter of the lesion, φ, is defined by the full-width half-maximum (FWHM) of the Gaussian sphere, i.e. φ = 2(2ln2)-1/2 σ. When a tumor prior Ct is included in the model, to ensure that the volume fractions add to unity (Ca + Cf + Ct = 1) at all locations r, we scale both Ca and Cf by (1 - Ct). The optical properties in the region of the tumor is then defined bywhere μ m0 is the baseline malignant tumor property (listed in Table 2) determined similarly to μ adi; parameter γ controls the simulated tumor contrast. We set γ to a range of values (including negative values) to validate our algorithm to the variation of tumor optical contrast. Note that μs’ values resulted from Eq. (3) are converted to scattering amplitude (Sa) and scattering power (Sp) for a multi-spectral reconstruction [3]. The developed digital breast phantom and associated data can be downloaded at http://openjd.sf.net/digibreast.2.2 Generation of the simulated optical measurements
The optical source and detector locations used the same spatial distributions as in our 2nd-generation combined optical/mammography imaging system [37], which has 96 continuous-wave (CW) sources (48 at 690 nm, 48 at 830 nm), 24 radio-frequency (RF) sources (for both wavelengths), 32 CW and 20 RF detectors. The source and detector locations are illustrated in Figs. 1(d) and 1(e), respectively.
A diffusion-equation based forward model [14] was used to generate simulated RF and CW optical measurements at 690 and 830 nm under various tumor settings. The diffusion equation was numerically solved on the forward mesh using our in-house finite-element solver – Redbird [38]. For each forward simulation, we added two pseudo-random noise sources to the model output to simulate realistic measurement noise:
where Φ0(s,d) is the fluence (mm−2) calculated by the diffusion model for the given source (s) and detector (d) pair; n 1 = A 1 × U × simulates the shot-noise and n 2 = max(|Φ0|) × A 2 × U simulates the electronic noise; U is a random variable in the standard normal distribution. The maximum amplitudes of n 1 and n 2, i.e. A 1 and A 2, were determined by the worst-case scenario of our instrument as reported in [37]. In this study, we set A 1 = 10−5 mm−1 (~0.1% noise for the highest signal and ~30% noise for the lowest signal) and A 2 = 10−5 (~100dB dynamic range).2.3 Improved structural prior-guided image reconstruction algorithm
The structural-prior guided reconstruction algorithm used in this study is described in [30]. Briefly, we first decompose the x-ray mammogram into tissue compositional (or volume fraction) maps, using a “fuzzy” segmentation approach. We then use these tissue compositional maps to build a regularization matrix, L, and run a Tikhonov-regularized Gauss-Newton reconstruction. For the k-th iteration of the Gauss-Newton reconstruction, we solve the following update equation (for under-determined cases [30]):
where µ = {HbO, HbR, Sa, Sp} is the solution vector defined over the mesh nodes; Θ is a block-diagonal matrix containing four (L T L) blocks along the diagonal; Jk = [J HbO, J HbR, JSa, JSp] is the Jacobian matrix at the k-th iteration; λ is the regularization parameter; I is the identity matrix; Φ denotes the measurement data; A denotes the forward operator.In this study, we developed an improved formula compared to [30] to enhance the computational efficiency. We first pre-compute the QR factorization of L, L = QR, and store R −1 in an upper-triangular matrix. For each Gauss-Newton iteration, we right-multiply the pre-computed R −1 to each sub-block in Jk to get Z = JR −1. Then the matrix multiplication ZZ T provides one of the diagonal blocks of Jk Θ−1 Jk T in Eq. (5) because
The pre-multiplication term Θ−1 Jk T in (5) is also simplified to R −1 Z T.This new formulation involves mostly multiplications by a triangular matrix, thus requiring less computation. Using this optimized formulation shortened both the reconstruction run-time and the peak memory requirement by a factor of 2 to 3.
In addition to our previous work where only two tissue compositions were used, in this study, we also used three tissue compositions – adipose, fibroglandular and tumorous tissues. Because our algorithm is general, for arbitrary numbers of tissue compositions, no additional change was needed aside from determining Ca, Cf and Ct as described earlier.
2.4 Validation of the structural-prior guided reconstructions
2.4.1 Fuzzy segmentation algorithms
We have implemented three segmentation approaches that can convert an x-ray mammogram into Ca and Cf maps. As reported in [30], the first approach (the “manual” approach) is to estimate the “pure” adipose and fibroglandular x-ray intensities by manually selecting regions-of-interest on the x-ray image, and then applying a linear mapping (Eq. (3) in [30]). In the second approach (the “dual-Gaussian” approach), we assume a Gaussian-mixture model [39] and fit the mammogram pixel intensity (χ) histogram using two Gaussian distributions, G a(χ) and Gf(χ), and then derive Cf using their cumulative distribution functions (CDF) by
In the third approach (the “threshold” approach), we first calculate the CDF of the histogram of the original x-ray mammogram [40]. Gray scales corresponding to a specified cutoff and (1-cutoff) values on the CDF curve are assigned as “pure” adipose and fibroglandular x-ray intensities [30], respectively. We then apply the same linear mapping as the Eq. (3) in [30] to convert any intermediate grayscale into fractional values between 0 and 1. Two cutoff values were tested, 0.2% and 2%. In this study, we evaluated the comparative performance of these segmentation methods in terms of 1) accuracy: comparing the segmentation errors against the reference, i.e. the fibroglandular volume fractions measured by the MicroDose SI system, and 2) effectiveness: comparing the errors of the reconstructed optical images against the simulation ground-truth.2.4.2 Regularization parameter settings
The regularization parameter, λ in Eq. (5), controls the strength of the priors and outcome image contrast [30]. Here we logarithmically increase λ from 10−6 to 10 with a step size of 101/4, and compare the recovered total hemoglobin concentration, HbT = HbO + HbR, images with the “ground-truth” defined in the breast phantom. The root-mean-square errors (RMSEs) are calculated for all reconstructions under the set of λ values, and the desired regularization is defined as the λ at which minimum RMSE is obtained. The same procedure is repeated for non-prior guided, two-composition-prior (using only Ca and Cf) guided, and three-composition-prior (using Ca, Cf and Ct) guided joint image reconstruction regimes with both uniform and locally refined meshes.
2.4.3 Simulations of tumors
A range of tumor contrasts (γ is set to ± 0.5, ± 1, 2, 3 and 4) and sizes (φ is between 4 mm to 2 cm) are simulated and subsequently reconstructed. After the image reconstruction, we extract the optical properties along a line profile, parallel to the x-axis (Fig. 1(c)), passing through the tumor center, and estimate the Weber contrast (R), defined by
where µ Ω represents the recovered optical property within 3 cm from the assumed tumor center, and µ b is the optical property in the surrounding background. Here we use the ground-truth (without lesion), i.e. µ(r) defined in Eq. (1), around the tumor regions as µ b.2.4.4 Impact of mesh resolution to tumor recovery
To assess the impact of the mesh resolution on the minimum tumor size that can be recovered, we rerun the reconstructions for all tumor sizes and contrasts at both assumed tumor region (Ωa and Ωf) using the locally refined meshes. The recovered tumor contrasts at all sizes are computed and compared against the results using the uniform mesh pair.
2.4.5 Influence of inaccurate prior information
It is important to test the susceptibility of our algorithm to faulty tumor priors. To this end, we intentionally define the tumor prior at the alternative site, i.e. for the tumor located in Ωa, we use a tumor prior centered at Ωf in the reconstruction, and vice versa. The reconstructed images using the discordant and concordant tumor priors are compared.
3. Results
3.1 Comparisons between fuzzy segmentation algorithms
The fibroglandular tissue volume fraction maps (Cf) derived from various fuzzy segmentation methods are compared against the MicroDose SI spectral measurement (ground-truth) in Fig. 2 . The RMSE between each segmentation method and the ground-truth fibroglandular volume fraction, as shown in Fig. 1(b), is marked in each sub-figure. The comparisons between the histograms of Cf in the compositional space using different segmentation approaches are compared in Fig. 2(e). Using the fibroglandular-adipose tissue priors (Ca and Cf) generated from each segmentation method, we ran two-composition-prior guided reconstructions and the RMSEs in HbT over 10 iterations are plotted in Fig. 3(a) . As a reference, we also plot the RMSEs using the measured Cf (i.e. ground-truth) as the prior. Recovered HbT image slices with and without the priors are compared against the HbT ground-truth in Figs. 3(b)-3(d).
Fig. 2 The comparison in accuracy of various fuzzy segmentation algorithms. Here we show the estimated fibroglandular volume fraction (Cf) using (a) the manual approach, (b) the dual-Gaussian approach, (c) the threshold approach with a 0.2% cutoff, and (d) the threshold approach with a 2% cutoff. The RMSEs between the estimated Cf and the ground-truth (Fig. 1(b)) are also superimposed on the figures in (a-d). The histograms of the estimated Cf using different approaches (color coded) are compared against the ground-truth (black) in (e).
Fig. 3 (a) Comparison in effectiveness of various fuzzy segmentation algorithms: (a) RMSE of HbT (μM) vs. reconstruction iterations; (b) HbT ground-truth (μM); (c) HbT (μM) recovered using two-composition-prior derived from the dual-Gaussian approach; (d) HbT (μM) recovered without using priors. In (b-d), black solid lines indicate Cf = 30% based on the ground-truth. All HbT images are extracted from the middle z-slice.
Instead of using the best-performing segmentation method, we pick the worst-performing one, i.e. the dual-Gaussian approach, for all subsequent simulations. This way, we can estimate the algorithm performance in the worst-case scenario and minimize the “inverse crime” by using different tissue compositional maps for the forward and inverse problems.
3.2 Selection of regularization parameters
In Fig. 4 , we show the RMSE vs. λ plots for three sets of reconstructions: 1) without structural priors, 2) with two-composition-prior and 3) with three-composition-prior, respectively, using the dual-Gaussian segmentation algorithm for a φ = 10 mm and γ = 2 tumor inside an adipose tissue background (Ωa).
Fig. 4 RMSE curves for reconstructions without tissue compositional priors (solid dots), and with two- (stars) or three- (circles) compositional priors using both uniform (blue) and locally refined (red) meshes. Priors are derived using dual-Gaussian segmentation algorithm and the tumor has γ = 2 in contrast and φ = 10 mm in size and is located in the adipose background.
3.3 Recovery of tumors
Using the identified λ value (0.032) based on the results in the above subsection, we ran a series of two- and three-composition-prior reconstructions using the simulated data for tumors over a range of sizes (φ) and contrasts (γ) located within either adipose or fibroglandular tissues. We plot a representative set (size φ = 4, 8, 12, 16, 20 mm, contrast γ = 1) of the recovered HbT line-profiles across the tumor centroids in Fig. 5 . The recovered tumor HbT contrasts (R) for all simulated lesion sizes and 4 representative contrast levels (γ = −1, 0.5, 2 and 4) are summarized in Figs. 6(a) and 6(b). In Figs. 6(c) and 6(d), the recovered HbT images of a small lesion (φ = 5 mm, γ = 2) surrounded by adipose tissue using two- and three-composition priors are plotted for comparison.
Fig. 5 The extracted HbT line-profiles across the (a-c) adipose-surrounded tumor and (d-f) fibroglandular-surrounded tumor, including (a,d) the ground-truth, and the profiles resulted from (b,e) two- and (c,f) three-composition-prior reconstructions. The blue dotted lines mark the boundaries of the averaging interval, i.e. Ω in Eq. (8), used in the contrast calculation.
Fig. 6 The recovered tumor contrasts (R) using the uniform mesh set for both (a) adipose- and (b) fibroglandular-surrounded tumors at various contrasts and sizes. Both three- and two-composition-priors are tested. Insets show the zoom-in views. Dash-dotted lines of the same color represent the true simulated contrast. Sample recovered HbT images (µM) of a small lesion (φ = 5 mm, γ = 2, adipose-surrounded) using two- and three-composition-priors are plotted in (c) and (d). Black solid lines mark the Cf = 30% contour on the ground-truth. Tumor is indicated by a dashed circle but with a radius twice the true size for easier visualization.
3.4 Mesh resolution and its impact on tumor reconstructions
We reran the above tests using the locally refined mesh set, and the recovered tumor contrasts are plotted in Figs. 7(a) and 7(b). The recovered HbT image for the above lesion is shown in (c).
Fig. 7 Comparisons between the uniform and refined meshes in the recovery of HbT contrasts (R) for tumors located in either (a) adipose or (b) fibroglandular vicinities. Insets show the zoom-in views. Markings are similar to those in Fig. 6. (c) The recovered HbT image (µM) of the same adipose-surrounded tumor as in Fig. 6(d), except that the refined mesh set was used.
3.5 Tolerance to faulty tumor priors
In Fig. 8 , we show the recovered HbT images reconstructed using both true and false tumor priors. In Figs. 8(a)-8(e), a 10 mm tumor of contrast γ = 2 is located in Ωa, while in Figs. 8(f)-8(j), the tumor is located inside Ωf. The reconstructed HbT images with the concordant and discordant priors are shown in (c,h) and (d,i), respectively. True tumor regions are marked by black dashed circles, whereas falsely assumed tumor priors are marked by white dotted circles.
Fig. 8 Test of tumor recovery with different tissue background and inaccurate prior information: (a-e) tumor is located in the adipose tissue and (f-j) tumor is in the fibroglandular tissue. In both cases, we show the (a,f) HbT ground-truth (µM), (b,g) compositional priors of the fibroglandular tissue, (c,h) recovered HbT images (µM) using concordant tumor priors, and (d,i) those using discordant tumor priors (white dotted circle). In all cases, the tumor has φ = 10 mm and γ = 2. In comparison, the two-composition-prior guided HbT images are shown in (e,j).
4. Discussions and conclusions
Based on Fig. 2, the “threshold” approach with a 0.2% cutoff showed the best segmentation accuracy among the four tested segmentation methods. Judging from the histograms in Fig. 2(e), the segmentation errors mostly arise from the overestimation of fibroglandular tissue content. Moreover, for linear segmentation algorithms such as the “threshold” and “manual” approaches, assuming the existence of “pure” tissue constituents results in a moderate “pooling” effect towards both ends of the histogram, which is absent from that of the ground-truth. Particularly, for the “threshold” method, the higher the cutoff value, the larger the portion of “pure” tissue constituents, as indicated in Fig. 2(e).
Comparing the results in both Figs. 2 and 3(a), we found that the accuracy of the fuzzy segmentation is directly related to how effective it is in terms of recovering the optical images when used as the structural prior – the errors of the recovered HbT decreased with decreasing segmentation errors. In other words, the more accurate the compositional priors, the higher the accuracy of the recovered optical images. Despite this difference, the RMSE variations in the recovered HbT images due to different priors were actually quite small, on the order of only ~0.02 µM on average (Fig. 3(a)). Comparing to the nearly two-fold increase in RMSE for the reconstructions without using a prior, all prior-guided reconstructions seem to produce consistently better optical images, especially in regions with poor optical fiber coverage, as shown in Figs. 3(b)-3(d). From a processing perspective, automated segmentation algorithms, such as the “threshold” and “dual-Gaussian” approaches, do not involve manual steps and can be adapted to the x-ray histogram characteristics that may vary between different mammographic systems, and are, therefore, more desirable.
From Fig. 4, the impact of λ on the prior and non-prior guided reconstructions of this particular phantom can be assessed. A trade-off can be observed here: when λ is below 10−4, the diminishing regularization results in rapidly increasing optical property estimation errors due to the contributions from the magnified measurement noise and modeling errors; on the other hand, larger λ values (>1) also result in increasing RMSE errors, likely a result of losing flexibility to fit the data as the prior starts to dominant the solution. Based on the results in Fig. 4, a balance was found roughly between 0.01 to 0.1, with no significant difference with regard to the meshes, i.e. numbers of unknowns when solving Eq. (5), nor to the priors used. For most tested reconstruction settings, a λ value around 0.032 appears to minimize the HbT estimation errors, and thus was selected for all our subsequent reconstructions.
A number of relevant findings can be obtained from the images in Fig. 6. First, using only fibroglandular and adipose tissue structures, our two-composition-prior algorithm was able, although with limited accuracy, to recover the tumor contrast – the recovery was more accurate in larger lesions (φ > ~16 mm) and in adipose rather than fibroglandular background. Second, the use of tumor priors (Ct) resulted in a 1.5- to 4-fold enhancement in the recovered tumor contrast by comparing the R values marked by circles with those by stars, indicating the importance of specifying potential tumor locations. Third, the reconstructions for an adipose-surrounded tumor appear to be more effective than for a fibroglandular-surrounded tumor, although the three-composition-prior guided reconstructions were less impacted by the tumor-to-background contrast. Furthermore, the recovered tumor contrasts in most cases underestimated the true contrast (dashed lines in Fig. 6); the smaller the tumor size, the more severe the underestimation. This is apparently a result of regularization and has been reported in previous studies [27]. Despite the underestimation, when the simulated optical tumor contrast is sufficient, e.g. γ ≥ 2, using three-composition-prior guided reconstructions, tumors as small as 5 mm were visually discernable from surrounding adipose tissues, as shown in Fig. 6(d). Finally, it appears that there was an effective lower bound on the minimum lesion size that could be recovered, and its value depended on both optical tumor contrast γ and the background tissue type. By visually inspecting the recovered images, a 10% or more deviation from the background (i.e. R ≥ 0.1) in a localized region appears to be discernable (as the case in Fig. 6(d)). Applying this criterion just for comparison purposes, the detection limit for three-composition-prior guided reconstructions of an adipose-surrounded tumor was approximately 10 mm for | γ | < 1, and 4 - 5 mm for | γ | ≥ 2; using the two-composition priors required roughly twice the tumor contrast to reach the same detection limit. If the tumor was located in the fibroglandular tissue, the limit increased to 8 mm for | γ | ≥ 1.
We did not include the plots for non-prior guided reconstructions in Fig. 6, but they closely tracked the results from the two-composition-prior guided reconstructions in almost all tested tumor sizes/contrasts. We believe this is because both methods do not contain tumor priors, and thus behave similarly in tumor recovery. Nonetheless, based on Figs. 3(a), 3(c) and 3(d), the use of normal tissue structure priors does significantly reduce the overall HbT error and improve the spatial details over the entire normal tissue regions.
By using the locally refined meshes, as shown in Fig. 7(a), the recovered tumor contrasts for small tumors (φ ≤ 8 mm) showed consistent improvement for adipose-surrounded tumors. This is further validated by inspecting the recovered HbT image shown in Fig. 7(c) of the same tumor case as in Figs. 6(c) and 6(d). However, such enhancement was not seen for fibroglandular-surrounded tumors. Moreover, the recovered tumor contrast using the refined meshes never overshot the ground-truth, as it did when using the uniform meshes, especially for tumors larger than 14 mm, indicating a more reliable and robust reconstruction. While using the refined mesh could lead to more accurate contrast estimates, only mild improvement is observed in lowering the tumor detection size limit: 8 mm for | γ | < 1. Moreover, smaller tumors residing within fibroglandular tissue still remain challenging with no improvement observed compared to results of reconstructions using the uniform mesh, as shown in Fig. 7(b).
The images in Fig. 8 further suggest the robustness of our algorithm. When an accurate tumor prior was provided, a significant increase in optical contrast was observed; however, when the tumor prior deviated from the truth, the algorithm was relatively insensitive to the false tumor priors and automatically fell back to the two-composition-prior guided reconstructions (Figs. 8(e) and 8(j)). As a result, no significant HbT contrast was observed near the falsely assumed tumor location.
In addition to the above tests using the worst-performing segmentation algorithm, i.e. the “dual-Gaussian” approach, we also ran the same analysis using the best-performing algorithm, i.e. the “threshold” approach with 0.2% cutoff. As we anticipated based on the results in Fig. 3(a), the recovered tumor contrast plots (not shown) followed closely to those presented in Figs. 6 and 7. The insensitivity to segmentation methods is a demonstration of robustness of our prior-guided reconstruction algorithm. We believe this is also one of the key reasons for its low susceptibility to false priors.
In summary, we have systematically tested our structural prior guided DOT reconstruction algorithm using a clinically derived complex digital breast phantom based on direct x-ray spectral measurements. We showed that the optical property estimation error can be reduced by almost 50% by using the structural priors, and the more accurate the structural priors, the more accurate the reconstructed optical images. By running reconstructions for a range of simulated tumor sizes and contrasts, we further demonstrated the importance of using tissue structural priors, especially tumor priors, in the reconstructions. For high-contrast (| γ | ≥ 2) adipose-surrounded tumors, combining the prior knowledge of the tumor location with the structural prior is shown to reduce the minimum detectable tumor size from 10 mm to 4-5 mm; for the fibroglandular-surrounded tumor, while this approach has also resulted in a 1.5 to 4 fold increase in the recovered tumor contrast, the minimum detectable tumor size was not changed. Moreover, we also found that using a refined reconstruction mesh could help recover smaller tumors residing within adipose tissue with improved image contrast, thus facilitating a more confident diagnosis. Finally, we demonstrated the robustness of our algorithm towards incorrect tumor priors. Only a concordant tumor prior can produce a significant increase in optical contrast. The distinct responses to concordant and discordant tumor priors allow our algorithm to potentially pinpoint the nature of “suspicious” regions when used during the image reading process. If automated, it can even serve as a computer aided detection, or CAD, method to automatically find suspicious tumors.
In the next steps, we will further fine-tune our segmentation algorithms to better approximate the measured tissue compositions. We will retroactively apply the algorithm with optimized parameters to the analysis of our clinical data set. Over the past 10 years, we have collected over 450 clinical measurements using our combined DOT and x-ray tomosynthesis imaging system [9,37], including over 200 patients with breast lesions (including malignant, solid benign, and cystic lesions). We expect this optimized reconstruction approach will lead to more accurate quantification of breast tumors. In addition, we will also expand this analysis to breasts of different density categories. We will specifically optimize our algorithm to help finding malignant tumors in the dense breasts, which is currently a key clinical challenge.
Acknowledgments
The authors acknowledge the Cooperative Research Matching Grand funding support, grand #2011D000334, from the Massachusetts Life Sciences Center, Massachusetts General Hospital Executive Committee on Research Interim Support Funding, and the NIH funding support under grants R01-CA97305, R01-CA142575, R01-GM114365 and 1S10RR023043.
References and links
1. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001). [CrossRef]
2. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]
3. R. Choe and A. G. Yodh, “Chapter 18. Diffuse optical tomography of the breast,” in Emergenging Technology in Breast Imaging and Mammography, (American Scientific Publishers, 2008) pp. 317–342.
4. V. Venugopal and X. Intes, “Recent advances in optical mammography,” Curr. Med. Im. 8(3), 244–259 (2012). [CrossRef]
5. P. Taroni, A. Torricelli, L. Spinelli, A. Pifferi, F. Arpaia, G. Danesini, and R. Cubeddu, “Time-resolved optical mammography between 637 and 985 nm: clinical study on the detection and identification of breast lesions,” Phys. Med. Biol. 50(11), 2469–2488 (2005). [CrossRef] [PubMed]
6. L. C. Enfield, A. P. Gibson, N. L. Everdell, D. T. Delpy, M. Schweiger, S. R. Arridge, C. Richardson, M. Keshtgar, M. Douek, and J. C. Hebden, “Three-dimensional time-resolved optical mammography of the uncompressed breast,” Appl. Opt. 46(17), 3628–3638 (2007). [CrossRef] [PubMed]
7. R. Choe, S. D. Konecky, A. Corlu, K. Lee, T. Durduran, D. R. Busch, S. Pathak, B. J. Czerniecki, J. Tchou, D. L. Fraker, A. Demichele, B. Chance, S. R. Arridge, M. Schweiger, J. P. Culver, M. D. Schnall, M. E. Putt, M. A. Rosen, and A. G. Yodh, “Differentiation of benign and malignant breast tumors by in-vivo three-dimensional parallel-plate diffuse optical tomography,” J. Biomed. Opt. 14(2), 024020 (2009). [CrossRef] [PubMed]
8. Q. Zhu, P. U. Hegde, A. Ricci Jr, M. Kane, E. B. Cronin, Y. Ardeshirpour, C. Xu, A. Aguirre, S. H. Kurtzman, P. J. Deckers, and S. H. Tannenbaum, “Early-stage invasive breast cancers: potential role of optical tomography with US localization in assisting diagnosis,” Radiology 256(2), 367–378 (2010). [CrossRef] [PubMed]
9. Q. Fang, J. Selb, S. A. Carp, G. Boverman, E. L. Miller, D. H. Brooks, R. H. Moore, D. B. Kopans, and D. A. Boas, “Combined optical and X-ray tomosynthesis breast imaging,” Radiology 258(1), 89–97 (2011). [CrossRef] [PubMed]
10. S. P. Poplack, T. D. Tosteson, W. A. Wells, B. W. Pogue, P. M. Meaney, A. Hartov, C. A. Kogel, S. K. Soho, J. J. Gibson, and K. D. Paulsen, “Electromagnetic breast imaging: Results of a pilot study in women with abnormal mammograms,” Radiology 243(2), 350–359 (2007). [CrossRef] [PubMed]
11. A. Cerussi, D. Hsiang, N. Shah, R. Mehta, A. Durkin, J. Butler, and B. J. Tromberg, “Predicting response to breast cancer neoadjuvant chemotherapy using diffuse optical spectroscopy,” Proc. Natl. Acad. Sci. U.S.A. 104(10), 4014–4019 (2007). [CrossRef] [PubMed]
12. R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: A case study with comparison to MRI,” Med. Phys. 32(4), 1128–1139 (2005). [CrossRef] [PubMed]
13. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]
14. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42(5), 841–853 (1997). [CrossRef] [PubMed]
15. A. H. Gandjbakhche, R. Nossal, and R. F. Bonner, “Resolution limits for optical transillumination of abnormalities deeply embedded in tissues,” Med. Phys. 21(2), 185–191 (1999). [CrossRef] [PubMed]
16. B. R. White and J. P. Culver, “Quantitative evaluation of high-density diffuse optical tomography: in vivo resolution and mapping performance,” J. Biomed. Opt. 15(2), 026006 (2010). [PubMed]
17. S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16(7), 5048–5060 (2008). [CrossRef] [PubMed]
18. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, and B. J. Tromberg, “Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain,” Opt. Lett. 30(11), 1354–1356 (2005). [CrossRef] [PubMed]
19. N. Ducros, C. D’Andrea, A. Bassi, G. Valentini, and S. Arridge, “A virtual source pattern method for fluorescence tomography with structured light,” Phys. Med. Biol. 57(12), 3811–3832 (2012). [CrossRef] [PubMed]
20. B. Pogue, S. Davis, F. Leblond, M. Mastanduno, H. Dehghani, and K. Paulsen, “Implicit and explicit prior information in near-infrared spectral imaging: accuracy, quantification and diagnostic value,” Philos. Trans. A Math, Phys. Eng. Sci. 369(1955), 4531–4557 (2011). [CrossRef]
21. V. Ntziachristos, A. G. Yodh, M. D. Schnall, and B. Chance, “MRI-guided diffuse optical spectroscopy of malignant and benign breast lesions,” Neoplasia 4(4), 347–354 (2002). [CrossRef] [PubMed]
22. B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, T. D. Tosteson, J. Weaver, S. P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures by using hybrid MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. U.S.A. 103(23), 8828–8833 (2006). [CrossRef] [PubMed]
23. P. C. Pearlman, A. Adams, S. G. Elias, W. P. Th. M. Mali, M. A. Viergever, and J. P. Pluim, “Mono- and multimodal registration of optical breast images,” J. Biomed. Opt. 17(8), 080901 (2012). [CrossRef] [PubMed]
24. B. Brooksby, H. Dehghani, B. Pogue, and K. Paulsen, “Near-infrared tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron. 9(2), 199–209 (2003). [CrossRef]
25. A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42(25), 5181–5190 (2003). [CrossRef] [PubMed]
26. B. Brooksby, S. Jiang, H. Dehghani, B. W. Pogue, K. D. Paulsen, J. Weaver, C. Kogel, and S. P. Poplack, “Combining near-infrared tomography resonance imaging to study in vivo and magnetic breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt. 10(5), 051504 (2005). [CrossRef] [PubMed]
27. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, C. M. Carpenter, S. Jiang, and K. D. Paulsen, “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15(13), 8043–8058 (2007). [CrossRef] [PubMed]
28. X. Intes, C. Maloux, M. Guven, B. Yazici, and B. Chance, “Diffuse optical tomography with physiological and spatial a priori constraints,” Phys. Med. Biol. 49(12), N155–N163 (2004). [CrossRef] [PubMed]
29. D. Hyde, E. L. Miller, D. H. Brooks, and V. Ntziachristos, “Data specific spatially varying regularization for multimodal fluorescence molecular tomography,” IEEE Trans. Med. Imaging 29(2), 365–374 (2010). [CrossRef] [PubMed]
30. Q. Fang, R. H. Moore, D. B. Kopans, and D. A. Boas, “Compositional-prior-guided image Reconstruction Algorithm for Multi-modality Imaging,” Biomed. Opt. Express 1(1), 223–235 (2010). [CrossRef] [PubMed]
31. E. Fredenberg, M. Lundqvist, B. Cederström, M. Åslund, and M. Danielsson, “Energy resolution of a photon-counting silicon strip detector,” Nucl. Instrum. Methods 613(1), 156–162 (2010). [CrossRef]
32. H. Ding and S. Molloi, “Quantification of breast density with spectral mammography based on a scanned multi-slit photon-counting detector: a feasibility study,” Phys. Med. Biol. 57(15), 4719–4738 (2012). [CrossRef] [PubMed]
33. E. Fredenberg, D. R. Dance, P. Willsher, E. Moa, M. von Tiedemann, K. C. Young, and M. G. Wallis, “Measurement of breast-tissue x-ray attenuation by spectral mammography: first results on cyst fluid,” Phys. Med. Biol. 58(24), 8609–8620 (2013). [CrossRef] [PubMed]
34. P. C. Johns and M. J. Yaffe, “X-ray characterisation of normal and neoplastic breast tissues,” Phys. Med. Biol. 32(6), 675–695 (1987). [CrossRef] [PubMed]
35. Q. Fang and D. Boas, “Tetrahedral mesh generation from volumetric binary and gray-scale images,” Proc. of IEEE Int. Symp on Biomed. Imaging (ISBI’09), 1142–45 (2009), URL: http://iso2mesh.sf.net
36. R. Bade, H. Haase, and B. Preim, “Comparison of fundamental mesh smoothing algorithms for medical surface models,” In Proc. of Simulation and Visualization. 289–304 (2006).
37. B. Zimmermann, M. Martino, A. Sajjadi, Q. Fang, D. Boas, and S. Carp, “A novel tomographic optical breast imaging system to simultaneously co-register x-ray tomosynthesis,” OSA BIOMED conference, (2014)
38. Q. Fang, S. A. Carp, J. Selb, R. Moore, D. B. Kopans, E. L. Miller, D. H. Brooks, and D. A. Boas, “A multi-modality image reconstruction platform for diffuse optical tomography,” OSA BIOMED conference (2008)
39. M. A. Balafar, “Gaussian mixture model based segmentation methods for brain MRI images,” Artif. Intell. Rev. 41(3), 429–439 (2014). [CrossRef]
40. J. C. Russ, The Image Processing Handbook, 6th ed. (CRC Press, 2011)
