2.1.

Spatial Frequency Domain Imaging System

A compact SFDI system, purchased from Modulating Imaging Inc., was used to quantitatively image optical parameters at four NIR wavelengths (658, 730, 850, and 970 nm) using spatially varying light. In series, 30 spatial frequencies uniformly distributed from $0{\mathrm{mm}}^{-1}$ up to $0.33{\mathrm{mm}}^{-1}$ were projected per wavelength onto each sample using high power light-emitting diodes (LEDs), a projection system and a digital micro-mirror device. Total data acquisition (DAQ) time was 10 min on average. Structured light patterns were generated using C# code and an ALP applications program interface, MI Acquire (v1.3.12), provided by Modulated Imaging Inc. (Irvine, California). A 12-bit CCD-based camera, coregistered with the projector, captured narrow bands of spectral reflectance using interference filters. Illuminating at a small angle to normal incidence and cross-linear polarizers placed at the source and detector minimized specular reflections. Figure 1 shows the fully integrated projection and camera subsystems mounted on a $z$-axis post.^1,37

Fig. 1

Compact SFDI mounted on a $z$-axis post in light-shielded container; imaging a calibration standard in reflectance geometry.

2.2.

Data Acquisition and Calibration

A planar and harmonically varying source illuminated each sample in reflectance geometry at normal incidence according to

Here, $q_0(z)$ represents the decay of the source as a function of depth in the media, $\alpha$ represents its phase offset, and $f_x=k_x/2\pi$ describes the spatial frequency of modulation in the $x$-direction. Assuming a linear medium and symmetry considerations, each sinusoidal source gave rise to a reflected intensity ($I_{\mathrm{AC}}$) with the same frequency and phase:¹

A three-phase demodulation was then employed to extract the frequency dependent modulation amplitude, $M_{\mathrm{AC}}$, from the reflected intensity and to remove ambient noise common to all three images:³⁸

The modulation amplitude of a siloxane titanium dioxide (${\mathrm{TiO}}_2$) reflectance standard with expected optical properties, $M_{\mathrm{AC},\mathrm{ref}}$, was measured at each imaging session for phantom-based calibration and to correct for spatial nonuniformity in the illumination and imaging systems. The reference standard was characterized by diffuse optical spectroscopy (DOS), which combines frequency domain photon migration with broadband spectroscopy.³⁹ Diffuse reflectance of the standard ($R_{d,\mathrm{ref}}$) was predicted from its expected optical properties using a forward model of light transport. Absolute reflectance of the sample ($R$) was subsequently calculated according to

Each spatial sampling point was calibrated and analyzed separately.

2.3.

Comparison of Analytical Models of Reflectance

Optical parameters were fit by minimizing the residual sum of squares between the measured modulation amplitude and its analytical solution in a semi-infinite, turbid medium, modeled according to diffusion¹⁶ and scaled-Monte Carlo simulations of light transport.²⁵ The former approximation assumes isotropic scattering and the latter, a Henyey-Greenstein phase function.¹ Rescaling a single, spatially and temporally resolved Monte Carlo simulation significantly reduced the computational burden associated with iterative estimation of optical parameters. However, computing its Fourier transformation was still time intensive in a minimization scheme. Consequently, scaled-Monte Carlo simulations of the modulation amplitude were stored in a 30-frequency, 4-wavelength, look up table (LUT) to reduce the computational burden associated with iterative estimation of spectral parameters. The LUT modulation amplitude was computed for absorption coefficients varying from 0 to $0.0999{\mathrm{mm}}^{-1}$ in steps of $1.9057E\text{-}04{\mathrm{mm}}^{-1}$ and scattering coefficients varying from 0.2495 to $2.1495{\mathrm{mm}}^{-1}$ in steps of $0.0019{\mathrm{mm}}^{-1}$. A LUT was not implemented for the diffusion approximation because computation of its forward analytical solution was direct; however, other groups have implemented this to improve speed of optical parameter rendering.¹ Cuccia et al. compared the accuracy of the diffusion approximation and scaled-Monte Carlo simulation of light transport using computer simulations for a range of optical pathlengths, and showed that for low spatial frequencies in a medium with an albedo approaching 1, the diffusion error remains less than 16%.¹ Here, model-data goodness of fit was evaluated in modulated reflectance measurements acquired from 47 patients in a HIPAA-compliant, prospective study, approved by Institutional Review Board for the protection of human subjects common to evaluate the applicability of these models to light scattering in surgical breast tissues. The adjusted $R^2$ coefficient of determination was used to quantify data variance explained by the model and the lowest adjusted Akaike information criterion (AIC) was used to rank model performance because it is based on the negative log likelihood of the residuals and consequently, provided a relative assessment of nonlinear model fit (the $R^2$ coefficient of determination is known to inadequately assess nonlinear regression because the total sum of squares is not equal to the regression sum of squares, as in the linear case).^40,41 A spectrally constrained optimization was explored with the diffusion approximation given a priori knowledge of the chromophore extinction spectra and the wavelength-dependence of the reduced scattering coefficient. However, an unconstrained optimization was observed to be more robust given the limited number of wavelengths sampled. The ability of recovered spectral parameters to distinguish breast tissue types in this large tissue population are evaluated elsewhere; here the optimization of model-data fit is presented.

2.4.

Optical Parameter Recovery

Spectral parameters were subsequently fit according to absorption by the endogenous tissue chromophores, oxygenated-hemoglobin (${\mathrm{HbO}}_2$), deoxygenated-hemoglobin (Hb), and water, and the wavelength-dependence of light scattering. Concentrations of ${\mathrm{HbO}}_2$ and Hb were recast as total hemoglobin ($\mathrm{HbT}={\mathrm{HbO}}_2+\mathrm{Hb}$) and percent oxygen ($\%{\mathrm{O}}_2={\mathrm{HbO}}_2/\mathrm{HbT}$). A power law dependence on wavelength was used to describe the Mie scattering spectrum: the scattering amplitude (A) and scattering power (b) were calculated relative to the reduced scattering coefficient of 1% intralipid at ${\lambda }_0=800\mathrm{nm}$⁴²

Adipose was not included in the basis set of fit chromophore extinction spectra because the four wavelengths sampled were insufficient to uniquely separate water and lipid absorption; an assumption about lipid absorption was necessary during optimization to prevent saturation of percent lipid or percent water at its bounds, as previously discussed by Mazhar.⁴³ Suspected lesions were primarily fibroglandular and pathology-assessed lesions were on average, 3% adipose, as determined prospectively by immunohistochemistry.

2.5.

Tissue Simulating Phantoms

A homogeneous, polymer-based phantom with inclusion wells varying in diameter (1.25 to 15 mm) was custom designed and constructed by INO Inc. (Biomimics phantom, Quebec City, Quebec, California) for contrast-detail analysis. The absorption and reduced scattering coefficients of the phantom at 800 nm were 0.008 and $1.006{\mathrm{mm}}^{-1}$, respectively; the refractive index was 1.52 and carbon black was used as an absorber. The inclusions were filled with dilutions of Intralipid, as shown in Fig. 2, to determine the minimum size of detectable scattering contrast in reduced scattering maps. Actual scattering contrast varied from 1% to 66% and was compared to measured scattering contrast; plotted along the center $x$-axis of each phantom according to:

Fig. 2

Contrast-detail curves for scattering heterogeneities that vary in size ($x$-axis) and contrast ($y$-axis); measured at (a) 658 nm, (b) 730 nm, (c) 850 nm, and (d) 970 nm.

A series of homogenous liquid phantoms with expected optical properties were imaged to test the absolute quantification of absorption and reduced scattering coefficients by SFDI. Phantoms were constructed from Intralipid-20% (Fresenius Kabi, Uppsala, Sweden), a fat emulsion used to mimic scattering in biological tissues, diluted with phosphate buffered solution (PBS) and porcine blood.⁴² PBS was used to dilute the blood to maintain blood cells at their physiological pH of 7.4. Scattering concentrations varied from 0.5% to 2.0% and the reduced scattering coefficient of each dilution was computed per wavelength according to Michels et al.⁴² Initial dilutions of porcine blood were prepared using measured hematocrit values assayed with a HemoCue© meter (Hb 201 system, HemoCue Inc., Cypress, California) for each sample, and the final solution values were serially diluted from 60 to 3.75 µM, the expected physiological range in breast tissue.⁴⁶ Expected optical properties were subsequently calculated using the forward model of light transport and assuming the chromophore extinction spectra.⁴⁴ The Pearson’s correlation coefficient, and coefficient of variation for all expected and recovered optical parameters were computed for model comparison. The coefficient of variation is just the standard deviation normalized to the distribution mean, and approaches 1 when two variables change together exponentially.

2.6.

Spatial Frequency Domain Imaging of Lumpectomies

Four fully intact and un-inked lumpectomy specimens were imaged immediately upon excision at the time of primary surgery to test the clinical translation of this new technology. The surgeon participated in each imaging session and placed sutures in the excised tissues to safeguard knowledge of its orientation until surgical inks were applied and the tissue was sent to pathology. Fifteen spatial frequencies uniformly distributed between 0 and $0.33{\mathrm{mm}}^{-1}$ were projected per wavelength onto each margin in under 2 min and spectral parameters were reconstructed using the scaled-Monte Carlo LUT. Tissues were extremely malleable, so only one to three margins were imaged per specimen.

3.1.

Contrast-Detail Analysis

Maps of the reduced scattering coefficient are shown in Fig. 2 for the contrast-detail phantoms and contrast line scans are plotted in Fig. 3(a) through 3(d). Detectable scattering contrast decayed asymptotically with inclusion size, as shown in Fig. 4(a): a 1.25 mm diameter inclusion was detectable when scattering contrast was at least 28%; the lowest detectable scattering contrast was 7% when the inclusion diameter was at least 10 mm. Recovered scattering contrast was observed to be lower than the expected scattering contrast, indicated by the dashed lines in Fig. 3(a) through 3(d), likely due to assumptions made about the optical properties of the inclusions. Figure 4(a) shows that absorption contrast, due to variations in phantom water content, was not detectable for inclusions less than 5 mm in diameter and absorption maps were characterized by a low signal-to-noise ratio (SNR). Figure 5(a) shows the reduced scattering coefficient maps for three phantoms with scattering inclusions placed at depths of 0 to 3 mm. True scattering contrast is mapped in Fig. 5(b); at least 33% scattering contrast was detectable at a 1.5 mm depth. The minimum detectable scattering contrast as a function of inclusion depth is plotted in Fig. 4(b) for the gelatin phantom; similar levels of detectable scattering contrast were observed in the polymer-based phantom. No absorption contrast was introduced into the gelatin phantom and absorption maps showed low signal to noise ratio (SNR) and minimal cross talk with scattering signatures, particularly with increasing inclusion depth.

Fig. 3

(a)–(d) Detected contrast in reduced scattering maps for each contrast detail phantom per wavelength (indicated by color); scattering contrast is plotted along the center $x$-axis of each phantom. Expected scattering contrast per wavelength indicated by dashed line.

Fig. 4

(a) Minimum detectable scattering contrast plotted as a function of inclusion diameter for the reduced scattering coefficient and absorption coefficient. (b) Minimum detectable scattering contrast plotted as a function of inclusion depth for the reduced scattering coefficient. Absorption contrast was not detectable for inclusion diameters less than 5 mm in these phantoms.

Fig. 5

(a) Reduced scattering maps (${\mathrm{mm}}^{-1}$) as a function of wavelength and depth for a gelatin phantom with true scattering contrast illustrated in (b) and homogenous absorption ($0.005{\mathrm{mm}}^{-1}$).

3.2.

Absolute Quantification of Optical Parameters

Recovered and expected optical parameters for Intralipid phantoms, serially diluted with porcine blood and extracted from the scaled-Monte Carlo LUT, are plotted in Fig. 6(a) and 6(d). The reduced scattering coefficient was systematically underestimated (on average, by 26%), but relative trends in the spectral scattering amplitude and slope were conserved. Systematic error in the absolute quantification of scattering was likely influence by the assumed optical parameters of the calibration phantoms, which for Intralipid, was derived by Michels⁴² using Mie theory and a different acquisition geometry. Spectral parameters extracted from the NIR optical coefficients are plotted as a function of Intralipid and hemoglobin concentration in Fig. 6(b) and 6(c) as well as in 6(e) and 6(f). Recovered hemoglobin values demonstrated insensitivity to scattering variations except at low concentrations and were within $+/-12\%$ of their expected values. The absorption coefficient was overestimated at visible wavelengths for low hemoglobin concentrations ($<15\mathrm{uM}$) because of the shorter pathlengths sampled. Water quantification was independent of hemoglobin changes and oxygenation levels approached 100% for all liquid phantoms (data not shown). Systematic under-estimation of absorption parameters at 970 nm was observed, likely due to the assumptions about the expected water absorption in the phantoms. Errors in Fourier Sampling theory may have also contributed minimally to amplitude underestimation, and consequently absorption overestimation, at low spatial frequencies.¹ The scattering coefficient demonstrated insensitivity to hemoglobin variations, except at the measured concentration limits (0 and 60 μM). The scattering amplitude increased linearly with scattering concentration (low coefficient of variation with Intralipid concentration: 0.26), but the scattering slope depended on both the particle size and number density (high coefficient of variation with Intralipid concentration: 1.36). All recovered absorption and scattering parameters showed strong correlation with their expected values ($>0.97$).

Fig. 6

(a) and (d) True (blue) and recovered (red) absorption and reduced scattering coefficients extracted from homogeneous, tissue simulating phantoms using scaled-MC simulations as a forward model. Intralipid concentrations varied from 0.5% 1.0% and 2.0%. Hemoglobin concentrations varied from 0 to 60 uM in six serial dilutions. Parameters recovered at all four wavelengths are shown in sequence. (b) and –(c) Spectral absorption parameters mapped as a function of hemoglobin concentration at three scattering concentrations (indicated by line color). (e) and (f) Spectral scattering parameters mapped as a function of Intralipid concentration for six hemoglobin concentrations (indicated by line color). Error bars indicate the standard deviation observed per $40\times 40$ pixel window.

3.3.

Forward Model Selection

Table 1 compares the diffusion approximation and scaled-Monte Carlo simulations of light transport according to model-data goodness of fit in a large cohort of surgical breast tissues. Three variations of the diffusion approximation were explored: a spectrally unconstrained model (a) and two spectrally constrained models (b) and (c), whose optimized parameters are listed in Table 4.2, to investigate the dependence of GOF on spectral constraints. The reduced scattering and absorption coefficients were constrained, respectively, to a power law dependence on wavelength and exponential extinction according to tabulated chromophores. Model (c) implements the most severe spectral constraints by linking the total concentration of hemoglobin to oxygen saturation (bounded between 0% and 100%) during optimization. GOF decreased with increasing spectral constraints, likely due to sparse spectral sampling. The $R^2$ coefficient of determination shows that the spectrally unconstrained diffusion model and scaled-Monte Carlo model performed comparably with regards to model-data goodness of fit, however the lowest adjusted-AIC suggests that the scaled-Monte Carlo model best minimized information loss in this population of surgical tissues. The correlation between the recovered and expected absorption and reduced scattering coefficients was greater than 0.97 for all models; however their coefficient of variation (CV) indicated some exponential trends in absorption estimates with absorber concentration for both the diffusion ($\mathrm{CV}=0.86$) and scaled-Monte Carlo models ($\mathrm{CV}=0.71$); likely because the high spatial frequencies probed are less sensitive to the longer photon pathlengths associated with absorption. Sparse spectral sampling may have also limited the accuracy of chromophore quantification. Recovered scattering coefficients varied linearly with their expected values for both models ($\mathrm{CV}=0.27$ to 0.28). The scaled-Monte Carlo model was selected because its recovered absorption parameters varied more linearly with hemoglobin concentration in liquid phantoms within a physiologic hemoglobin range.

Table 1

Goodness of fit evaluated in breast tissue for three variations of the diffusion approximation and a scaled-Monte Carlo forward model, according to their average adjusted R2 coefficient and the adjusted Akaike information criterion (AIC).

3.4.

Challenges Associated with the Clinical Translation of SFDI

Spectroscopic images of four surgical margins are shown in Fig. 7: column A shows two wire-localized fibroadenomas (0.7 to 0.8 mm in diameter) in a fibrocystic background and column B shows four small fibroadenomas (0.2 to 0.5 cm in diameter) in an otherwise complete pathologic response to neoadjuvent chemotherapy. Column C shows a nearly complete pathologic response to chemotherapy, except one remaining high-grade tumor focus (0.2 cm) in the fibrocystic background. The lymphadenectomy associated with this partial responder is shown in column D; 11 out of 15 lymph nodes were positive for cancer. Oxygenation was the most apparent gauge of tumor response to therapy in these four lumpectomy specimens; oxygen levels in the lymph nodes filled with tumor were $\sim 40\%$ lower than those observed in pathologically responsive tumors and fibroadenomas. Variation in optical parameters with time delay after excision reported by Bydlon suggest that oxygen saturation levels of hemoglobin vary nonlinearly within 30 min of resection; whereas, the reduced scattering coefficient and total hemoglobin concentration showed little change over time.⁴⁷ Consequently, oxygenation could be a valuable diagnostic parameter if imaged in the surgical cavity or immediately at the time of resection. High spectroscopic scattering was observed from fibrocystic tissues (columns A and C) and a lower scattering slope distinguished two large fibroadenomas within a fibrocystic background, as shown in column A. Water content was elevated and oxygenation levels were lower in the upper left quadrant of the partial responder (column C), suggesting that this quadrant was the site of the high-grade tumor focus. The primary challenges observed with the clinical translation of SFDI were edge artifacts, which presented as elevated absorption at the tissue edge, a typical location of tissue profile decline, and orienting spectroscopic images with pathology to validate margin status.

Fig. 7

Photograph and spectral parameter maps sampled from four lumpectomy margins: (a) fibroadenomas ($n=2$) with diameters of 0.7 to 0.8 mm surrounded by fibrocystic disease, (b) complete pathologic response to neoadjuvent chemotherapy and two small fibroadenomas (diameters of 0.2 to 0.5 mm); (c) partial response to neoadjuvent chemotherapy with residual high grade tumor focus (0.2 mm) and fibrocystic disease; (d) lymphadenectomy associated with case (c). Cancer is present in 11/15 nodes. Skin and muscle in the image field appear white and deep red, respectively in the tissue photograph.