2.1.

Sd-SFDI Process

An overview of the sd-SFDI process can be seen in Fig. 2(a). To image the samples, we used a custom built SFDI system, diagramed in Fig. 2(b). A DLP Lightcrafter Evaluation Module digital micromirror device with built-in LEDs (Texas Instruments, Dallas, TX) was used to project sinusoidal light patterns of varying spatial frequencies and phases onto a sample. Based on experimental testing, a minimum of 12 pixels on the DMD was necessary to recreate accurate grayscale sinusoidal patterns. Due to the limited resolution of the DMD, a plano-convex singlet lens ($f=15\mathrm{cm}$) was placed in front of the DMD’s stock projection lens to decrease the size of the projected pixels, enabling the system to achieve spatial frequencies exceeding $1.3{\mathrm{mm}}^{-1}$. The system would otherwise be limited to only $\sim 0.4{\mathrm{mm}}^{-1}$ based on the minimum projection distance using the included projection optics. The sinusoidal patterns were projected onto the sample at an angle orthogonal to the patterns’ direction, and all samples were flat or flattened when imaged. This setup helped reduce specular reflections⁴⁵ and made measurements less sensitive to variations in the samples’ height.⁴⁶ By mitigating changes in the projected patterns’ spatial frequency as a function of distance from the projection axis, this orthogonal setup also made measurements less sensitive to the incidence angle and variations thereof over the length and width of the sample.⁴⁶

Fig. 2

Sub-diffuse spatial frequency domain imaging method. (a) The full flow diagram for imaging a tissue sample and rendering a heatmap of its optical properties. This example uses four spatial frequencies. The tissue is imaged with the three-phase measurement method at various frequencies; the images are demodulated and calibrated for each frequency; and the resulting calibrated image cube is input into a machine learning model which calculates the optical properties. (b) To image the tissue, sinusoidal light patterns of varying spatial frequencies are projected onto the tissue, and the reflected light is measured with a camera. The axis of projection is orthogonal to the direction of the sine waves to decrease specular reflections and sensitivity to sample height variations. (c) A neural network similar to the example shown here is used to transform the reflectance spectra of each pixel into an optical property measurement. Each node’s inputs are aggregated in a weighted sum and passed through an activation function to yield the node’s output.

A 5-megapixel monochrome CCD camera (piA2400-17gm Basler, Ahrensburg, Germany) was used to collect images of the light reflected from the sample.⁴⁷ We employed the standard three-phase measurement method, by which measurements were acquired at three phases for each spatial frequency. These measurements were then demodulated.¹¹ This process was repeated for 26 spatial frequencies (0.01, $0.025{\mathrm{mm}}^{-1}$, and 24 spatial frequencies evenly spaced between 0.05 and $1.25{\mathrm{mm}}^{-1}$) and 3 different wavelengths (450, 530, and 620 nm), resulting in an uncalibrated [$2050\times 2048\text{pixels}\times 26$ spatial frequencies] image cube per wavelength, each pixel of each image cube containing a unique demodulated intensity spectrum, $M(f)$. A calibration phantom was measured with this same method to find $M_{\mathrm{ref}}(f)$ per pixel per wavelength. Each pixel of each uncalibrated image cube was then calibrated using the equation $R(f)=M(f)/M_{\mathrm{ref}}(f)*R_{\mathrm{ref}}(f)$,¹¹ resulting in a calibrated image cube for each wavelength.

The spatial frequencies used in this study include 6 frequencies in the diffuse domain ($f\le 0.2{\mathrm{mm}}^{-1}$), 5 frequencies in the intermediate domain ($0.2{\mathrm{mm}}^{-1}<f<0.5{\mathrm{mm}}^{-1}$), and 15 spatial frequencies in the sub-diffuse domain ($f\ge 0.5{\mathrm{mm}}^{-1}$). The spatial frequencies in the diffuse domain were included as they improved estimation of ${\mu }_s^{\prime }$ (Supplement 1 in the Supplementary Material), which affects $R(f)$ in both the sub-diffuse and diffuse domain.² This also meant that optical property estimation had to account for possible variations in ${\mu }_a$, which affects $R(f)$ in the diffuse domain.¹¹ The upper bound of $1.25{\mathrm{mm}}^{-1}$ was chosen because higher frequencies showed a poor signal-to-noise ratio due to limitations in the dynamic range and sensitivity of the camera.⁴⁸ Optionally, additional demodulation processing¹¹ can be applied to the images to include the DC frequency, $f=0{\mathrm{mm}}^{-1}$, but validation tests showed no increase in prediction quality when including the DC frequency (data not shown), likely due to the frequencies very close to $0{\mathrm{mm}}^{-1}$ which are already included. As such, we did not include images at $f=0{\mathrm{mm}}^{-1}$ in our procedure.

An ANN similar to that developed by Naglič et al. was developed to predict optical properties from sd-SFDI reflectance spectra [Fig. 2(c)]. The model was a multi-layer perceptron with two hidden layers and an output layer,³⁷ implemented with Keras and configured to work on OpenCL.³⁶ The model, which has several weights that must be trained by fitting to observed data, ultimately behaves as an empirically trained analytical equation which maps the $R(f)$ and $g_1$ of a pixel to an optical property value. As such, it computes very quickly. Two instantiations of this model were created and trained for ${\mu }_s^{\prime }$ and $\gamma$, respectively, allowing each model to more finely tune to its respective optical property (memory was not a considerable limitation). While one model with two outputs may have worked instead, using one model constrains the trained weights of the intermediate layers to have values that lead to accurate results for two outputs simultaneously, and only the final layer would be able to contain weights unique to each property. Using two distinct models rids us of this constraint. These models were used to render ${\mu }_s^{\prime }$ and $\gamma$ heatmaps at all three wavelengths.

2.2.

Generating Simulated Dataset and Training the Model

A simulated dataset was generated to train the ANN and the benchmark method’s semi-empirical equation. A diffuse Monte Carlo program developed by Hennessy et al.⁴⁹ and Alerstam et al.⁵⁰ was modified to include the parameter $\gamma$. This was achieved by developing a novel phase function sampling method. In this method, the phase function’s inverse CDF equation used for inverse transform sampling is replaced with an inverse CDF look-up table using methods similar to those introduced by Naglič et al.³² In order to allow for simulation across a wide range of the phase function parameters $\gamma$ and $g_1$, our look-up table is generated with the modified Henyey-Greenstein²³ (MHG) function when $\gamma$ and $g_1$ fall within MHG’s valid range ($\gamma \le 1+g_1$³¹) and is otherwise generated with the GK.⁵¹ We approximated GK’s valid range to be $\gamma >1+0.6g_1$ and $\gamma <3^{g_1}$³⁸ (Supplement 1 in the Supplementary Material), so our method cannot sample phase functions with $\gamma \ge 3^{g_1}$. This is not a concerning limitation given the range of phase function parameters seen in our phantoms and in biological tissue.^{2,21,26,32–35} Our method uses a novel mapping which relates the $\gamma$ and $g_1$ values to the GK parameters³⁸ (Supplement 1 in the Supplementary Material) to calculate the GK parameters when using the GK phase function. It should be noted that the GK phase function has an analytical inverse CDF⁵² that can be used in tandem with our mapping in place of the inverse CDF look-up table, but using look-up tables for both phase functions allowed for more seamless switching between GK and MHG, the latter of which does not have an analytical inverse CDF.³²

This hybrid MHG-GK phase function sampling method allows for simulations of $R(f)$ across a wide range of phase function parameters encompassing the $\gamma$ and $g_1$ values seen in both our phantoms and in biological tissue.^{2,21,26,32–35} The method’s ability to sample phase functions with a specified $\gamma$ and $g_1$ combination, enabled through our GK mapping and MHG’s analytical relation to $\gamma$ and $g_1$,²³ allows for the simulation of a reference phantom’s spectrum $R_{\mathrm{ref}}(f)$. This is necessary for calibrating measurements of physical data so that the calibrated data can be input into models trained on simulated data.¹¹

We verified that the phase function of the hybrid model was a good match for the phase function of the calibration phantom as calculated with Mie theory (Supplement 1B in the Supplementary Material). The modified Monte Carlo program with hybrid phase function sampling was then used to generate a training dataset, a validation dataset, and the calibration spectrum. Datasets were generated for combinations of the optical properties seen in Table 1, chosen to cover the range of optical properties seen in biological tissues^{2,21,26,32–34} as well as our tissue phantoms. A wide array of ${\mu }_a$ values was included to account for variations in ${\mu }_a$ values when estimating ${\mu }_s^{\prime }$. Table 1 also shows the optical properties of the phantoms whose measurements served as the experimental dataset for testing the model. The labels of the datasets, which the models would be trained to predict, were ${\mu }_s^{\prime }$ and $\gamma$. The features, which the models would use to make these predictions, comprised the reflectance values $R(f)$ at the 26 spatial frequencies as well as the value of $g_1$. Table 1 shows the total number of samples used in each dataset. Combinations that fell outside the upper bound of the GK range (i.e., $\gamma \ge 3^{g_1}$) were discarded during simulation, shrinking the total number of samples in the training and validation dataset.

Table 1

Parameters used for training, validation, and experimental datasets. Formats of “[x,y]n”, signify n elements in the range of [x,y].

The value of $g_1$ is assumed to be known ahead of time, allowing for its use as a feature. This is possible because $g_1$ is known for the phantoms and is assumed to be constant at 0.9^45,53 for biological tissues. More specifically, the $g_1$ of soft mammalian tissue falls within the range of 0.8 to 0.95,³⁵ and $R(f)$ remains relatively invariant to changes in $g_1$ at $g_1$ values $>0.8$.⁵⁴

We note that the training data did not entirely encompass the range of our testing data, which is an opportunity for future improvement with the model. However, the model still performed well when applied to testing data whose optical properties fell outside the range of the data the model was trained on.

The model was trained on the training dataset, and the validation dataset was used to determine the stopping point. More specifically, the model was set to train for 2000 epochs or until the mean squared error of the model on the validation dataset failed to improve for 100 epochs in a row.

While outside the scope of this study, we note that our sd-SFDI method could potentially also find ${\mu }_a$ values. The model already includes $R(f)$ values at diffuse spatial frequencies in its inputs and is trained on simulations with a variety of ${\mu }_a$ values. However, the Monte Carlo simulations that the ANN was trained on used a semi-infinite geometry, where $x$ and $y$ are infinite in space, and $z$ extends to negative infinity, which means that the ANN assumes the sample has infinite thickness. Although this assumption holds at sub-diffuse frequencies, where the penetration depth of the photons is shallow,²⁶ the assumption does not hold when considering our tissue samples at diffuse frequencies. At diffuse frequencies, the penetration depth is in some cases deeper than the thickness of our tissue samples,²⁵ which range from around 1 to 2 mm in thickness. Fortunately, the impact this has on our model’s ${\mu }_s^{\prime }$ calculations, which only partially relies on $R(f)$ values at diffuse spatial frequencies, appears to be minimal, as shown in Supplement 2 in the Supplementary Material. However, the prediction of ${\mu }_a$ is likely to be adversely affected by low sample thickness, so further research would be required to examine and account for this impact before attempting to predict ${\mu }_a$. Alternatively, a future study could use thicker tissue samples.

2.3.

Creating and Measuring Tissue-Mimicking Phantoms

Tissue-mimicking phantoms were created using polystyrene beads suspended in water (Polysciences Inc., Warrington, PA) with mean bead diameters of $d_1=0.0878\mu \mathrm{m}$ (SD: $0.01\mu \mathrm{m}$), $d_2=0.19\mu \mathrm{m}$ (SD: $0.01\mu \mathrm{m}$), and $d_3=0.99\mu \mathrm{m}$ (SD: $0.03\mu \mathrm{m}$). We prepared two dilutions of each bead distribution stock solution in distilled water: one to achieve a value of ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$ at 530 nm and the other to achieve a value of ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$ at 530 nm, as calculated using Mie theory.⁵⁵ This resulted in six bead phantoms which covered $\gamma$ values ranging from 0.97 to 2.17 and ${\mu }_s^{\prime }$ values ranging from 1.11 to $5.39{\mathrm{mm}}^{-1}$ to cover the ranges of values previously reported for various biological tissues.^{2,21,26,32–34} The optical properties of the phantoms at each wavelength were quantified using Mie Theory.⁵⁵ All phantoms had negligible absorption (assumed to be $0.001{\mathrm{mm}}^{-1}$), as they were made with non-absorbing beads and distilled water. The six phantoms were measured at three wavelengths, resulting in eighteen total sd-SFDI image cubes.^47,56 No cross polarization was used during sd-SFDI measurement as to avoid rejecting sub-diffuse photons.⁵⁷

One of these six phantoms was measured before the rest of the phantoms to serve as a calibration reference, thus creating the reference spectrum $M_{\mathrm{ref}}(f)$ for each pixel and wavelength. When simulating this phantom to find $R_{\mathrm{ref}}(f)$ at each wavelength, ${\mu }_a$ was set to $0.001{\mathrm{mm}}^{-1}$ to represent negligible absorption, as a value of 0 caused problems with the simulator.

2.4.

Measuring Skin Tissue Samples

Mohs surgery skin tissue samples were obtained in collaboration with the dermatologic surgery faculty at Dell Medical School. This study was approved and the informed consent requirement was waived by the Institutional Review Board at The University of Texas at Austin and the Seton Healthcare Family. We used four samples from three patients. Excision sites included the ear, the ankle, and the nose. The samples were frozen, the top layer of each sample was removed and processed using standard H&E staining methods, and digital microscopic images of the slides were captured and stitched.⁵⁸ The thicknesses of the removed layers ranged roughly from 0.5 to 3 mm. Regions of interest delineating epidermis, dermis, and basal cell carcinoma (BCC) were outlined by a board certified dermatologist and fellowship trained Mohs surgeon (MCF) from Dell Medical School. The remainder of the sample was thawed and imaged using our custom SFDI set-up, using the same method and reference phantom that was used for measuring the phantoms. During imaging, tissue samples were positioned on the stage such that they fell within the same area the reference phantom was positioned in during imaging, but outside the area of the specular reflection that the reference phantom had. The tissue samples were flattened under a 1-mm glass slide to mitigate additional specular reflections. Marked H&E images of the tissues were co-registered with the SFDI heatmaps and used as a mask for each tissue subtype so that the optical properties of their pixels could be compared.

It is important to note that the co-registration between digital histology images and optical property heatmaps is not exact, as is likewise reported in the previous studies.^2,6 The frozen sample from which the top layer is removed for histological analysis must be thawed for SFDI, resulting in slight shape changes to the sample. However, the heatmaps and H&E images are approximately the same shapes. The marked H&E images were manually rotated, scaled, and overlaid onto the heatmaps in MATLAB, and the markings were manually translated into the digital masks, keeping in mind other landmark features of the images to guide the making of the masks (Supplement 1C in the Supplementary Material).

Another limitation is that the histology image only captures the top-most layer, whereas the diffuse frequencies of the SFDI image are sensitive to deeper layers of the tissue.⁵⁶ Moreover, since the top layer of the sample is removed for histology before the sample is imaged with sd-SFDI, this top layer itself is not captured in the sd-SFDI images.

2.5.

Testing the Model

Average $R(f)$ spectra were extracted from the 15 phantom sd-SFDI image cubes by averaging the spectra from select pixels of the image cubes at each wavelength (Supplement 1B in the Supplementary Material). To quantify the accuracy of the trained ANN model, the model was applied to these 15 averaged $R(f)$ spectra along with their corresponding $g_1$ values to predict the optical properties of each phantom at each wavelength. The mean absolute relative error of the results for each optical property was calculated using the equation $|\overline{\epsilon }|=1/n*\sum _i^n(y_i-{\hat{y}}_i)/y_i$.

The model was then applied to the SFDI image cubes of the phantoms to render ${\mu }_s^{\prime }$ and $\gamma$ heatmaps of the phantoms. The SFDI image cubes of the skin tissue were then processed using the trained ANN to render heatmaps of ${\mu }_s^{\prime }$ and $\gamma$. The skin tissue was assumed to have an anisotropy factor of $g_1=0.9$.⁴⁵ Software timing tools were used to measure the frame rate at which the model was able to render these heatmaps.

In tandem, ${\mu }_s^{\prime }$ and $\gamma$ predictions were rendered using the semi-empirical non-linear fitting method first reported by Kanick et al.²⁶

3.1.

Model Speed and Accuracy on Phantoms

The results of the model when tested on the experimental phantom spectra and image cubes are shown in Fig. 3. Looking at select sub-diffuse optical property heatmaps [Figs. 3(a) and 3(b)], we see that the model is able to correctly predict the optical properties of these phantoms, whose true values are shown in white text beneath the samples. Moreover, the model is able to make these predictions over a wide field and render optical property heatmaps from the sd-SFDI image cubes. One limitation to note, though, is that the heatmaps contain small artifacts from specular reflections in the liquid phantoms. This is likely due in part to not being able to use cross polarization. We avoided using this technique because it could reject sub-diffuse photons,⁵⁷ but it is often used to filter out specular reflections.⁴⁵ In order to mitigate the degradation of the inverse model’s performance caused by these artifacts, the specular reflections were avoided during all quantitative analysis of phantom heatmaps. This was done by selecting a section from each phantom that was free from specular reflections (Supplement 1B in the Supplementary Material). In future experiments, specular reflections may be reduced further by subtracting measurements of water in a deep well with dark walls from both the sample and reference phantom measurements.^2,45

Fig. 3

Performance of the neural network across five homogeneous phantoms with known optical properties and three wavelengths. Image cubes comprising 26 sd-SFDI images were transformed into optical property heatmaps for each phantom at each wavelength using the trained deep learning model, resulting in 15 heatmaps. (a), (b) Select ${\mu }_s^{\prime }$ heatmaps and select $\gamma$ heatmaps, respectively, from the results, with the true value for the phantoms written in white underneath. Note these selections span multiple wavelengths and phantoms. The model accurately produces wide-field optical property heatmaps over the range of values tested. (c), (d) The accuracy of the model when run on an average spectrum from each of the phantom image cubes, where the $x$ axis is the true optical property, the $y$ axis is the predicted optical property, and the dotted green line represents unity. The mean absolute relative error was 6.8% for ${\mu }_s^{\prime }$ and 3.6% for $\gamma$.

Quantifying the accuracies further, Figs. 3(c) and 3(d) show the predicted values plotted against the true values, with the solid line representing unity. The mean absolute relative errors of the ${\mu }_s^{\prime }$ predictions and $\gamma$ predictions were 6.8% and 3.6%, respectively, which is within the range of error values previously reported in the literature.^2,6,26 When run on simulated versions of these same phantom spectra, the model performed with an accuracy of $\left|\overline{\epsilon }\right|=0.77\%$ for ${\mu }_s^{\prime }$ and $|\overline{\epsilon }|=0.94\%$ for $\gamma$, on par with the accuracy that Naglič et al.’s model achieved over a narrower range of $\gamma$ values.³⁶

Comparing the speed of this method to the non-linear fit method previously used for sd-SFDI (Table 2), the ANN is able to render optical property heatmaps from $2050\times 2048\text{pixel}$ image cubes in significantly less time. The frame rate of the non-linear fit method is only 0.0005 fps, i.e., a frame rendering time of 2000 s per frame. In contrast, the ANN run on the same device has a frame rate of 2.32 fps, i.e., a frame rendering time of 0.43 s per frame, which meets the requirement for real-time imaging.¹⁴ It is important to note that the ANN method is implemented with parallelization, allowing it to process multiple pixels simultaneously, while the non-linear fit method is not. In order to compare the speeds of the two methods independent of their parallelization, the speed of each process was examined when run on a single pixel. The non-linear fit method processed a single pixel with an average speed $0.013\mathrm{s}/\text{pixel}$, while the ANN accomplished this in only $0.00187\mathrm{s}/\text{pixel}$. Thus even without parallelization, the ANN still outperforms the non-linear fit method by an order of magnitude.

Table 2

Comparison of the frame rate for non-linear fit method versus artificial neural network.

3.2.

Tissue Type Comparison

Some highlighted results of the tissue experiments can be seen in Fig. 4. We focus on ${\mu }_s^{\prime }$ and $\gamma$ at 530 nm in this figure because these two properties showed the greatest quantitative differences between the tissue subtypes as determined by visual inspection of all quantitative figures. A comprehensive figure with all results can be seen in Supplement 3 in the Supplementary Material. Figure 4(a) shows zoomed in regions of interest from the tissue samples for each tissue subtype. Dashed lines circle the regions of interest on H&E images, ${\mu }_s^{\prime }$ heatmaps, and $\gamma$ heatmaps for each tissue type. Areas of contrast can be seen in the heatmaps that align with the boundaries of the tissue subtypes.

Fig. 4

(a) Selected regions of the H&E images, along with corresponding regions of the ${\mu }_s^{\prime }$ and $\gamma$ heatmaps at 530 nm, for each of the tissue subtypes. Areas of contrast which align with the marked regions of interest can be seen in the heatmaps. (b) Box and whisker plots of the optical properties of the pixels from each tissue subtype combined across all four tissue samples. BCC has lower ${\mu }_s^{\prime }$ than adipose and dermis, while adipose has higher $\gamma$ than dermis and BCC. (c) Two-dimensional scatter plot of the mean values of the optical property for each sample’s tissue subtypes ± one standard deviation. The tissue subtypes fall into three distinct clusters, with $\gamma$ separating adipose from BCC and ${\mu }_s^{\prime }$ separating Dermis from BCC.

Adipose tissue exhibits high ${\mu }_s^{\prime }$, likely due to connective tissue septa which are particularly dense in subcutaneous fat in areas of the face and head,⁴¹ from which 3 of our 4 samples were excised. Adipose tissue also exhibits high $\gamma$, indicating that the large adipocytes result in a high ratio of the adipose tissue’s largest particles to smallest particles.²⁴ Dermis exhibits a relatively high ${\mu }_s^{\prime }$ and a low $\gamma$, which may be caused by its collagen fibers and relatively small constituent particles.⁵⁹ BCC tends to have a lower ${\mu }_s^{\prime }$ than normal skin tissue subtypes, echoing trends seen in the previous research.^41,43,44 It also exhibits low $\gamma$ values, which may be associated with its high density of small nuclei.²

To examine the sub-diffuse optical property differences between tissue subtypes, Fig. 4(b) shows box and whisker plots of the combined pixels from the tissue subtypes across all four tissue samples. In these data, while ${\mu }_s^{\prime }$ values of BCC show a fair amount of overlap with those of adipose and dermis, the inner quartiles of ${\mu }_s^{\prime }$ values for BCC show almost no overlap with those of adipose and dermis. Meanwhile, $\gamma$ seems to separate adipose from dermis and BCC, as its inner quartiles show no overlap with those of dermis and BCC.

If these results prove to be consistent over a larger dataset, viewing them in the context of previous studies would suggest that $\gamma$ is a more reliable feature than ${\mu }_s^{\prime }$ for discriminating between adipose and BCC. Salomatina et al. conducted a study using spatially resolved diffuse measurements with wavelengths ranging from 370 to 1600 nm. The study found that adipose tissue in skin excised from the face, neck, and scalp exhibited trends similar to our results, with ${\mu }_s^{\prime }$ higher in adipose tissue than in BCC. However, this study also found that adipose tissue in skin excised from the back exhibited low ${\mu }_s^{\prime }$ values, on par with those of BCC, due to its thin connective tissue.⁴¹ McClatchy et al., in a study on breast tissue using sd-SFDI at wavelengths of 658, 730, and 850 nm, found that adipose tissue from breast also exhibits low ${\mu }_s^{\prime }$ values. However, this tissue still exhibited high $\gamma$ values similar to what is seen in our results.² These results indicate that $\gamma$ could be more consistently higher in BCC than in adipose tissue, whereas ${\mu }_s^{\prime }$ may not always be higher in BCC than in adipose tissue.

These findings suggest that using both of these properties in tandem could separate all three tissue subtypes from each other. To explore this idea further, Fig. 4(c) shows a two-dimensional scatter plot of the mean ${\mu }_s^{\prime }$ value versus the mean $\gamma$ value from each tissue’s adipose, dermis, and BCC regions, ± one standard deviation. The optical properties of each tissue subtype form small clusters with overlap within standard deviations but separated means.

Figure 5 examines the results of a select tissue sample. Looking at a single-tissue sample over a wide range is more representative of a clinical scenario, in which a medical professional may use sub-diffuse imaging to help discriminate between tissue subtypes. Figure 5(a) shows a reconstructed white light image of this tissue sample, created using the demodulated DC intensity images at the red, green, and blue wavelengths as channels in an RGB image. This is not a perfect reconstruction, as the individual wavelengths were not white-balanced, but it serves to approximate a white light photo of the tissue. Figure 5(b) shows the marked H&E image of this sample. Figures 5(c) and 5(d) show examples of the sd-SFDI reflectance images at a diffuse spatial frequency ($f=0.1{\mathrm{mm}}^{-1}$) and a sub-diffuse frequency ($f=0.7{\mathrm{mm}}^{-1}$), respectively, taken at a wavelength of 530 nm. The microstructural variety present within the tissue sample is not apparent when examining the white light image or the reflectance maps. However, looking at the sub-diffuse optical property heatmaps of the sample in Figs. 5(e) and 5(f), we see that the sample contains a wide range of optical property values indicating vast differences in microstructure. Moreover, looking at the tissue subtype masks that the marked histology image was used to form [Fig. 5(g)], we see that some of the contrast seen in the heatmaps aligns with the tissue subtype regions of interest. Finally, looking at a two-dimensional scatter plot of ${\mu }_s^{\prime }$ versus $\gamma$ for all pixels from these regions of interest [Fig. 5(h)], we see that the optical properties of the pixels from these subtypes fall into their own clusters, though these clusters do overlap.

Fig. 5

A representative tissue sample. (a) Reconstructed white light image of the sample. (b) Digital histology slide marked with regions of interest for the three skin tissue subtypes focused on for this study: adipose, dermis, and BCC. (c), (d) Demodulated and calibrated reflectance image for a diffuse and a sub-diffuse spatial frequency. Both images are taken at 530 nm. (e), (f) Optical property heatmaps of this sample for ${\mu }_s^{\prime }$ and $\gamma$, each at 530 nm. The images reveal structure not seen in the white light image. (g) Regions of interest from the marked histology slide translated onto the shape of the co-registered heatmaps. The regions of interest align with areas of contrast seen in the heatmap images. (h) ${\mu }_s^{\prime }$ versus $\gamma$ scatter plot of all the pixels from these regions of interest. Pixels from subtypes fall into optical property clusters, albeit with overlap.

In general, the heatmaps of the samples frequently showed contrast that corresponded roughly with the regions of interest, although this contrast did not go so far as to consistently create clear, obvious outlines of these regions of interest. The scatter plot shows clear trends between the optical properties of the tissue and the subtype of the tissue, supporting the hypothesis that sub-diffuse optical properties could be useful in identifying and discriminating between skin tissue subtypes.