2.1.

SFDI: Theory and Instrumentation

We previously developed a handheld SFDI device for breast imaging, described previously.³⁵ Briefly, the device consists of a monochrome CMOS camera (Texas instruments, Dallas, Texas) and digital light projector (LightCrafter, Texas instruments, Dallas, Texas) enclosed in acrylic housing. The original RGB LEDs of the projector were replaced with three high power surface-mount LEDs with peak emission wavelengths, $\lambda$, of 662, 735, and 859 nm (Luminus Devices, Sunnyvale, California). The two dichroic mirrors in the projector were replaced with ones with cutoff wavelength of 697 and 801 nm (FF697-SDi01 and FF801-Di02, Semrock optics, Rochester, New York). The original collimation lenses in front of the LEDs were also replaced (49874, Edmund Optics, Barrington, New Jersey). Crossed polarizers (LPVIS050 and LPNIRE100-B, Thorlabs, Newton, New Jersey) were added to reduce specular reflection. A glass optical window of diameter 50 mm (WG12012-B, Thorlabs, Newton, New Jersey) was used to make contact with the skin and four bar-style load cells were used to record the normal force applied to the tissue (Phidgets, Calgary, Alberta, Canada). Pattern projection and image acquisition from the camera was performed automatically by a custom MATLAB application implemented on a laptop computer. Real-time feedback of force applied was displayed on the computer screen in the form of a line graph. This was monitored by the operator to aid in applying a consistent pressure during the compression period.

Sinusoidal illumination at spatial frequency $f_x=0.1{\mathrm{mm}}^{-1}$ was projected at three equally spaced phase shifts for each wavelength, $\lambda$. The nine patterns (three phases, three wavelengths) were projected sequentially, with three blank frames included in the sequence to serve as a marker to identify the start of the sequence. Projection of the patterns was synchronized to the camera exposure with a frame rate of 40 fps. Each set of 12 frames, representing one multi-wavelength SFDI measurement, required a total of 300 ms for acquisition. One SFDI measurement was captured every 2 s. Following image acquisition, image demodulation was performed in real time. The zero frequency magnitude, $M(0)$, and the magnitude at $f_x=0.1{\mathrm{mm}}^{-1}$, $M(0.1)$, were calculated on a pixel-by-pixel basis according to Eqs. (1) and (2):²⁸

Calibration images were taken from a polydimethylsiloxane reference phantom of known optical properties, to obtain $M_{ref}(f_x)$. The reference phantom was assessed by a commercial frequency domain near-infrared spectroscopy system (OxiplexTS, ISS, Champaign, Illinois) and found to have ${\mu }_a$ of 0.047 and $0.052{\mathrm{cm}}^{-1}$ (at 690 and 830 nm, respectively) and ${\mu }_s^{\prime }$ of 7.1 and $5.3{\mathrm{cm}}^{-1}$. To estimate optical properties at the three $\lambda$ of interest (662, 735, and 859 nm), ${\mu }_{\mathrm{s}}^{\prime }(\lambda )$ was assumed to follow a power law decay function as provided in Eq. (3)²

2.2.

Human Subjects Imaging Protocol

Human subject testing was approved by the Carnegie Mellon Institutional Review Board under protocol STUDY2015_00000046. Thirteen healthy volunteers were recruited through the Pitt+Me online portal.³⁶ Six returned for a second imaging session, for a total of 19 imaging sessions. All subjects were female and ranged in age from 22 to 61 years with no history of diabetes or vascular disease. Enrolled subjects comprised a range of skin tones, as detailed below. Informed consent was obtained from all subjects at the beginning of each visit. A Fitzpatrick skin score was then assessed via a questionnaire adapted from Eilers et al.³⁷ The two questions were posed verbally to the subject and the experimenter would record her answers and assign a Fitzpatrick score based on the criteria described in Eilers et al.³⁷ The questions were “If after several months of not being in the sun, you stayed outdoors for about 1 hour at noon for the first time in the summer without sunscreen, what would happen to your skin? Would it become pink/red, irritated, tender, or itchy?” and “Over the next 7 days, would you develop a tan or notice your skin becoming darker?” The resulting Fitzpatrick scores of the study population are provided in Table 1. Fitzpatrick score of I corresponds to the lightest skin tone and VI to the darkest. Fitzpatrick score was unavailable for two subjects, which were estimated to be type III or lower.

Table 1

Distribution of Fitzpatrick skin scale scores for 11 of the 13 study participants. The remaining two participants are estimated to be type III or lower.

The experimental setup is shown in Fig. 1. The handheld device was used to apply local compression to the breast by pressing down against the ribcage at various locations on the breast. Subjects were lying in supine position while pressure was applied with the handheld device at two locations: one lateral to the nipple and one medial to the nipple. Five measurements were unusable due to technical failure of the imaging device, resulting in 71 measurements total. Each measurement consisted of 1 min of baseline measurement while the device was in contact with skin with minimal pressure, then one minute of sustained target pressure of 7.6 kPa, followed by a return to baseline pressure for one minute of recovery [Fig. 1(c)]. Petroleum jelly was used between the glass and the skin to ensure even contact and reduce trapped air. Over all measurements, an average pressure of $6.40\mathrm{kPa}\pm 0.85\mathrm{kPa}$ was reached during compression.

Fig. 1

(a) Subject positioning and direction of pressure applied with handheld imager. (b) Approximate location of two imaging locations on each breast. Image adapted from: Ref. 38 (c) average pressure applied over all 3-min measurements, shaded region indicates standard deviation ($n=71$).

2.3.

Signal Processing

2.4.

Processing of Dynamic Breast Compression Data

After optimizing the value of MI for each subject from their baseline measurement, the hemodynamic changes induced by compression of the breast tissue were calculated. This was achieved by calculating ${\mu }_{a,sub}(x,y,\lambda ,t)$ and ${\mu }_{s,sub}^{\prime }(x,y,\lambda ,t)$ on a pixel-by-pixel basis. Because of the presence of specular reflections and other artifacts in some images, pixels were eliminated if ${\mu }_{a,sub}(x,y)$ at any $\lambda$ or any $t$ differed by more than eight standard deviations from the average over a region $\sim 20\mathrm{mm}$ square in the center of the image. $\mathrm{HbO}(x,y,t)$ and $\mathrm{Hb}(x,y,t)$ were then computed from ${\mu }_{a,sub}(\mathrm{x},\mathrm{y},\lambda ,t)$ and converted to $\mathrm{tHb}(x,y,t)$ and ${\mathrm{StO}}_2(x,y,t)$ as described in Sec. 2.1. The circular field of view was partitioned into eight regions of interest (ROIs) and pixel values were averaged within each region to produce $\mathrm{tHb}(t)$ and ${\mathrm{StO}}_2(t)$ time traces. To assess stability in hemoglobin concentration changes during the baseline period and eliminate those regions with baseline drifts, linear fits were performed to the ${\mathrm{StO}}_2(t)$ trace for the time period $t=10$ to $t=55\mathrm{s}$. Baseline noise was assessed by calculating the root mean squared deviation of the data from the linear fit. ROIs with baseline slope of magnitude greater than 0.045 percentage points per second or variation $>0.4\%$ points were rejected. Linear fits were also performed for the time period $t=65$ to $t=113\mathrm{s}$ to obtain the rate of change of tissue oxygen saturation during compression.

3.1.

Monte Carlo Lookup Tables

To demonstrate the effect of the epidermis layer absorption on diffuse reflectance, four Monte Carlo LUTs at varying values of ${\mu }_{a,epi}$ are shown in Fig. 2. With the same isolines plotted for each, increasing ${\mu }_{a,epi}$ has the effect of shifting the table to lower values of $R_d(0)$ and to a lesser extent lower values in $R_d(0.1)$ while compressing the table to a smaller range of $R_d$ values in both dimensions. This represents decreased precision in quantifying optical properties at high ${\mu }_{a,epi}$, as small errors in the $R_d$ map lead to larger errors in ${\mu }_{a,sub}$ and ${\mu }_{s,sub}^{\prime }$. A given MI and $d_{epi}$ corresponds to a ${\mu }_{a,epi}$ value for each of the three wavelengths of interest. As determined empirically from in vivo data and the optimization procedure described in Sec. 2.3.3, the four LUTs represent a range from low to high melanin concentration in ${\mu }_{a,epi}(662)$.

Fig. 2

Example of two-layer LUTs with top layer absorption (${\mu }_{a,epi}$) fixed at four different values. For all, top-layer scattering is fixed at $16.89{\mathrm{cm}}^{-1}$ (corresponding to reduced scatting at 735 nm). Isolines for 14 increments of ${\mu }_{a,sub}$ are shown from 0.01 to $0.27{\mathrm{cm}}^{-1}$. Isolines for ${\mu }_{s,sub}^{\prime }$ are shown in 17 increments from 4 to $20{\mathrm{cm}}^{-1}$.

3.2.

Extraction of Melanin Index in Simulated Data

The forward simulation inputs are summarized in Fig. 3(a). With $d_{epi}$ of 0.011 cm, MI range from 0 to 0.16 corresponds to a ${\mu }_{a,epi}(662)$ range from 0 to $14.6{\mathrm{cm}}^{-1}$. For each MI, analysis with an LUT matching the ${\mu }_{a,epi}(\lambda )$ corresponding to the ground truth MI was found to reproduce ${\mu }_{a,sub}$ and ${\mu }_{s,sub}^{\prime }$ to within 1.4%. To estimate MI based on the optimization procedure in Sec. 2.3.3, we calculated optical properties using LUTs based on a range of MI values. This analysis showed that measurements of ${\mu }_{a,sub}(\lambda )$, and thus of tHb and ${\mathrm{StO}}_2$, were highly sensitive to the assumed value of MI. These trends are shown in Fig. 3(b) for two values of input MI (namely, 0.01 and 0.1.) The solid vertical lines indicate the actual MI value of the forward simulation. The dotted vertical lines correspond to the estimated MI based on minimum MSE obtained when fitting hemoglobin extinction coefficients to measured ${\mu }_{a,sub}$ for extraction of hemoglobin concentrations.

Fig. 3

(a) Illustration of the varying absorption properties of the eight multi-wavelength MC forward simulations, having $d_{epi}$ and ${\mu }_{s,sub}^{\prime }$ matched with the LUT. (b) Demonstration of the optimization of MI for simulated data, showing the effect of assumed MI on recovered tHb and ${\mathrm{StO}}_2$ values for two simulations of differing ground truth MI (0.01 and 0.1) and the MSE obtained from fitting the hemoglobin concentrations to ${\mu }_{a,sub}(\lambda )$. The MI value that produces the lowest MSE is marked with a dotted vertical line and the ground truth MI of the simulation is marked with a solid vertical line. (c) Percent error in tHb and ${\mathrm{StO}}_2$ quantification using the layered model (with MI obtained from lowest MSE) compared with a homogeneous model, showing drastically reduced error when the layered model is used.

These results illustrate that the point of minimum MSE slightly over-estimates the ground truth value of MI. This overestimation is reflected in the computation of underlying tHb and ${\mathrm{StO}}_2$, illustrated in Fig. 3(c). The mean error of ${\mathrm{StO}}_2$ extraction was found to be 2.9%, increasing from 1.6% at $\mathrm{MI}=0$ to 4.2% at $\mathrm{MI}=0.16$. The mean error for tHb was found to be $-1.9\%$, increasing from $-0.9\%$ at $\mathrm{MI}=0$ to $-3\%$ at $\mathrm{MI}=0.16$. Though error for the layered model is still dependent on MI, it is drastically reduced compared to the error from the homogeneous model, where the error in ${\mathrm{StO}}_2$ from $\mathrm{MI}=0$ to $\mathrm{MI}=0.16$ ranges from $-0.8\%$ to $-103.7\%$. and error in tHb ranges from 1% to 228.3%. Also plotted in Fig. 3(b) is the adjusted MSE, calculated as $\mathrm{MSE}/(1-{\mathrm{StO}}_2)$. This adjustment was developed based on the results of the in vivo data, detailed below. For simulated data, adjusted MSE yields the same estimate of MI as unadjusted MSE for all cases.

Error quantification for simulations where $d_{epi}$ and scattering differ from the properties assumed by the two-layer model are summarized in Fig. 4 for MI values up to 0.16. The upper limit of MI in this analysis was determined empirically, exceeding the greatest MI value calculated for the healthy subjects in Sec. 3.3. While a mismatch in scattering properties and layer thickness introduce erroneous results, all errors are lower than 11% and the majority of cases are lower than 5%. Human epidermal thickness is known to vary considerably, with one study reporting a range of $\sim 50$ to $100\mu \mathrm{m}$ for skin on breast (female), 90 to $165\mu \mathrm{m}$ for skin on abdomen and 100 to $175\mu \mathrm{m}$ for skin on the back,³⁹ and another reports a range of $\sim 50$ to $90\mu \mathrm{m}$ for forearm.⁴⁰ The thicknesses considered for this analysis (between $\sim 50$ and $170\mu \mathrm{m}$) were chosen to cover these ranges for epidermal thickness. These results indicate that though the two-layer model assumes a fixed epidermis thickness and scattering properties, this optimization procedure retains validity in estimating subcutaneous tHb and ${\mathrm{StO}}_2$ even as epidermis properties differ from their assumed values.

Fig. 4

Results of sensitivity analysis showing error in quantification of tHb and ${\mathrm{StO}}_2$ for mismatch in (a) epidermis thickness; (b) scattering amplitude; and (c) scattering power for four selected values of MI. For all rows, a negative mismatch represents a value which is lower than that assumed by the model, and a positive mismatch represents a value higher than that assumed by the model.

3.3.

Extraction of Melanin Index in Human Data

The effect of varying assumed MI on the goodness of fit of hemoglobin concentrations to ${\mu }_{a,sub}(\lambda )$ is examined for in vivo breast data. This trend, along with the effect on estimated tHb and ${\mathrm{StO}}_2$, is shown in Fig. 5 for two representative measurements from subjects of different Fitzpatrick scores.

Fig. 5

Example of the effect of assumed MI on recovered tHb and ${\mathrm{StO}}_2$ values for two subjects of different Fitzpatrick type. The third row shows the effect on MSE obtained from fitting the hemoglobin concentrations to ${\mu }_{a,sub}(\lambda )$, and the bottom row shows the adjusted MSE, obtained from $\mathrm{MSE}/(1-{\mathrm{StO}}_2)$. The MI value that produces the lowest adjusted MSE is marked with a dotted vertical line.

Similar to the trend observed in the simulated data, increasing assumed MI produced lower estimates of tHb and higher estimates of ${\mathrm{StO}}_2$. Unlike the simulated data, the relationship of MSE to MI was not unimodal but monotonically decreasing for all measurements, exhibiting the lowest value at the point corresponding to ${\mathrm{StO}}_2$ of 100%. Values above 100% or below 0% are considered undefined, as they would correspond to negative chromophore concentrations. This may be owing to the shape of the extinction spectrum of deoxyhemoglobin, which exhibits a dip around 740 nm. In all measurements, the measured value of ${\mu }_{a,sub}(735)$ is greater than that which would be predicted by the calculated hemoglobin concentrations alone. The failure of the measured ${\mu }_{a,sub}(\lambda )$ to conform to the shape predicted by the extinction spectra of hemoglobin may be due to the presence of other chromophores, such as water and lipid, which are not accounted for in this study. Regardless of its origin, the mismatch between observed and predicted ${\mu }_{a,sub}(735)$ leads to relatively large MSE when the fit yields physiologically reasonable values of Hb concentration, and much lower MSE for extremely low Hb concentrations, because the spectrum of HbO does not exhibit the dip in this range and a better fit is obtained.

Based on the apparent bias against high values of Hb, an adjusted MSE was introduced, in which MSE is divided by the corresponding percentage of Hb at that MI, ($1-{\mathrm{StO}}_2$). When this adjustment is applied, a minimum is present which corresponds to an average ${\mathrm{StO}}_2$ value over all measurements of $68\pm 11\%$ ($\text{mean}\pm \mathrm{std}$) and average tHb concentration of $19.2\pm 5\mu \mathrm{M}$. These ranges are physiologically reasonable based on literature for healthy breast.². The MI corresponding to minimum adjusted MSE is identified for each measurement. These values are averaged for each subject, and plotted against subject’s Fitzpatrick skin type in Fig. 6. All subjects of Fitzpatrick score I, II, and III received an optimal MI of 0, the single type IV subject received 0.003, and subjects of Fitzpatrick V and VI had a mean MI of 0.115 with a standard deviation of 0.052 (range: 0.014 to 0.166). Our results indicate that the algorithm is sensitive to higher MI values, as seen for participants of Fitzpatrick types V and VI exhibited higher MI, whereas the algorithm applied to participants of type I, II, III, and IV showed no sensitivity to MI, indicating that a homogenous reconstruction may be applicable.

Fig. 6

Relationship of Fitzpatrick skin score with MI, where MI is selected to produce the lowest error in fitting hemoglobin concentrations to ${\mu }_{a,sub}(\lambda )$. Each diamond represents one subject and error bars are standard deviation of the individual measurements. Due to the repeat imaging session by some subjects, $n$ ranges from 4 to 8 measurements. Types I and IV represent one subject each. Types III, V, and VI represent three subjects.

Using the adjusted MSE analysis, the optical properties could be extracted for all subjects. Specifically, each subject’s average estimated MI value (obtained from multiple measurements) is used to calculate baseline optical properties for all of that subject’s measurements. Baseline optical properties are also calculated for every measurement using the homogeneous LUT. Box plots of tHb, ${\mathrm{StO}}_2$, scatter amplitude $A$, and scatter power $b$ before and after correction are shown in Fig. 7 for two groupings of Fitzpatrick skin type: I, II, III, and IV (lighter skin tone) and V and VI (darker skin tone).

Fig. 7

Boxplots of optical properties obtained with homogeneous (top row) and layered (bottom row) inverse models for trials on Fitzpatrick I, II, III, and IV subjects (lighter skin tone) and subjects of Fitzpatrick V and VI (darker skin tone). The layered model result is based on the adjusted MSE described above and does not require a priori assumption of absorption coefficient of either layer. Differences between groups in the layered model values were assessed with two-sample $t$-tests, and a statistically significant difference is found to remain for tHb ($p=0.012$) and scatter power $b$ ($p<0.001$). For the homogeneous model values, $p<0.001$ was found for all four parameters.

Optical properties derived from the homogenous model reveal separation between groups not only for tHb and tissue saturation but also for scattering power and amplitude. Adjusted MSE-based correction brings the two groups closer together for each property considered. For the layered model results, the mean and standard deviation of tHb concentration was found to be $17.7\pm 4.6\mu \mathrm{M}$ for Fitzpatrick I, II, III, IV, and $21\pm 5\mu \mathrm{M}$ for Fitzpatrick V,VI, both in line with that reported by Grosenick et al.,² ranging from 12.6 to $19.2\mu \mathrm{M}$ tHb for healthy breast.² Similarly, using the layered model, ${\mathrm{StO}}_2$ for the two groups was $67.9\pm 9.3\%$ and $68\pm 13.6\%$, respectively, comparable to values reported by Grosenick et al.² of 67.7% to 74%. Layered model-derived scatter power, $b$, for the two groups was $0.90\pm 0.24$ and $0.67\pm 0.18$ compared to Grosenick et al. reporting values from 0.58 to 0.99. We also evaluated whether age is a contributing factor to the extracted optical data. Our data suggest (results not shown here) that there is a greater spread of optical properties with age when using the homogenous model. The greater difference in the homogenous assumption is likely explained due to higher average age of the Fitzpatrick V and VI subjects (48 years) than the Fitzpatrick I, II, III, and VI subjects (34 years). After correction with the layered model, no differences in optical properties were found with age.

3.4.

Breast Compression Dynamic Results

The tHb response to localized tissue compression with our handheld SFDI device was found to be consistent with tissue blanching and reperfusion, with tHb decreasing upon initiation of compression onset and recovering quickly after compression release. The response in ${\mathrm{StO}}_2$ was more variable and was therefore analyzed spatially according to eight ROIs as shown in Fig. 8. Out of 568 regions available for analysis, 132 of these were rejected due to lack of stable ${\mathrm{StO}}_2$ baseline ($n=70$) or for optical artifacts ($n=63$) such as specular reflection encompassing more than 60% of that region’s pixels. In the remaining 435 regions, a trend toward decreasing ${\mathrm{StO}}_2$ was observed during compression but with considerable variation. While the trend in ${\mathrm{StO}}_2$ sometimes differs within different regions of one measurement, inter-measurement and inter-subject variability was found to be even greater. Two representative measurements are shown in Fig. 8, one representing primarily decreasing ${\mathrm{StO}}_2$ trends in all eight regions and one representing primarily increasing ${\mathrm{StO}}_2$ trends. Plots of ${\mathrm{StO}}_2$ response to compression for all measurements are available in Figs. S1–S19 in the Supplementary Material. The red lines in Fig. 8 show a linear fit to the ${\mathrm{StO}}_2$ time trace during compression.

Fig. 8

Change in ${\mathrm{StO}}_2$ before during and after compression for two measurements taken from different subjects on different regions of the breast. The eight ROIs of each measurement are illustrated as gray shaded areas with an axis overlaid plotting the average change from baseline of ${\mathrm{StO}}_2$ in that region (black line), with a linear fit during the compression period overlaid in red. The missing axis represents a region rejected for unstable baseline or artifacts. Inset: $x$ and $y$ axis scale, which is consistent within each subject. The $y$ axis has units of $\mathrm{\Delta }{\mathrm{StO}}_2$ (%), and $x$ axis has units of seconds. Tick marks at 60 and 120 s mark the beginning and end of the compression period.

The overall variability of the slope of the ${\mathrm{StO}}_2$ responses (pp/s) across all ROIs is represented by the histogram in Fig. 9.

Fig. 9

Histogram of distribution of ${\mathrm{StO}}_2$ slopes during compression for 470 ROIs. ${\mathrm{StO}}_2$ slope is expressed in units of percentage points per second (pp/s). The median slope (marked with black vertical line) occurs at $-0.011\mathrm{pp}/\mathrm{s}$.

To illustrate the range of ${\mathrm{StO}}_2$ responses, the ROIs were segregated into three groups: those with a pronounced decreasing trend, those of intermediate slope, and those of pronounced increasing trend. The average responses of these three groups, shown in Fig. 10, indicate that the ${\mathrm{StO}}_2$ responses can be quite variable between different regions and across subjects, whereas the tHb response shows a more consistent decrease during compression. The group averages correspond to $n=145$, 145, and 145 ROI averaged time traces.

Fig. 10

Illustration of average change from baseline of (a) ${\mathrm{StO}}_2$ and (b) tHb for ROIs grouped according to ${\mathrm{StO}}_2$ slope during compression. Shaded region indicates standard deviation. From left to right, $n=145$, 145, and 145 ROIs. tHb exhibits a consistent U-shaped response across all three groups.

The ROIs were further analyzed in terms of the Fitzpatrick score of the subject they were obtained from. Again splitting the group into Fitzpatrick scale I, II, III, IV and V and VI, we found the ROIs trending down in ${\mathrm{StO}}_2$ contain 118 from I, II, III, IV, and 27 from V and VI, the ROIs exhibiting little change in ${\mathrm{StO}}_2$ contain 106 and 39, respectively, and the ROIs exhibiting increasing ${\mathrm{StO}}_2$ contain 66 and 79. While the distribution is not equal in all trends, there is not a clear distinction between trends and Fitzpatrick scale.