Imaging sub-diffuse optical properties of cancerous and normal skin tissue using machine learning-aided spatial frequency domain imaging

Abstract. Significance: Sub-diffuse optical properties may serve as useful cancer biomarkers, and wide-field heatmaps of these properties could aid physicians in identifying cancerous tissue. Sub-diffuse spatial frequency domain imaging (sd-SFDI) can reveal such wide-field maps, but the current time cost of experimentally validated methods for rendering these heatmaps precludes this technology from potential real-time applications. Aim: Our study renders heatmaps of sub-diffuse optical properties from experimental sd-SFDI images in real time and reports these properties for cancerous and normal skin tissue subtypes. Approach: A phase function sampling method was used to simulate sd-SFDI spectra over a wide range of optical properties. A machine learning model trained on these simulations and tested on tissue phantoms was used to render sub-diffuse optical property heatmaps from sd-SFDI images of cancerous and normal skin tissue. Results: The model accurately rendered heatmaps from experimental sd-SFDI images in real time. In addition, heatmaps of a small number of tissue samples are presented to inform hypotheses on sub-diffuse optical property differences across skin tissue subtypes. Conclusion: These results bring the overall process of sd-SFDI a fundamental step closer to real-time speeds and set a foundation for future real-time medical applications of sd-SFDI such as image guided surgery.

and the ratio of cell sizes. 1 This behavior can be characterized using optical properties such as the reduced scattering coefficient (μ 0 s ), absorption coefficient (μ a ), and the phase function backscatter parameter (γ). 2 These optical properties have been found to be useful in many medical applications such as detecting burn wound severity, 3 monitoring blood occlusions, 4,5 and aiding in cancer diagnostics. [6][7][8][9] Methods to measure such optical properties have been developed and improved. In fiber optic based approaches, photons are injected at one location, and the remitted photons are detected at short distances away and analyzed to determine the optical properties. 10 Although small fiber optics set-ups are effective, these point-based techniques have limited use in wide-field applications. To solve this issue, spatial frequency domain imaging (SFDI) has been developed, in which structured light is projected onto tissue and the light reflected is measured from the entire surface at once. 11 Using sinusoidal light patterns of various spatial frequencies/phases and demodulating the reflected intensity measured at each frequency, SFDI can produce a demodulated intensity spectrum MðfÞ for each pixel of the imaged tissue. Calibrating these spectra yields a reflectance spectrum RðfÞ for each pixel. Optical property heatmaps can then be rendered from this grid of reflectance spectra. [11][12][13] Recently, efforts have been made to advance the speed of SFDI up to real time, i.e., one frame per second, 14 which involves increasing the speed of both the measurement of SFDI images and the rendering of optical property heatmaps from these images. Many methods have been developed which can significantly reduce the acquisition time needed for capturing wide-field SFDI images. These include methods which capture images at multiple wavelengths simultaneously 15,16 as well as methods which reduce the number of SFDI snapshots needed to form an image. 14,17 There have likewise been many advancements in increasing the speed of processing these measurements and rendering optical property heatmaps. These methods include optimized look-up tables, 18 random forest regressors, 19 and deep learning networks. 4,[20][21][22] These advancements, while impressive, have only been implemented for diffuse scattering reflectance. As such, they cannot be used to find the phase function backscatter parameter γ, a parameter which has been shown to be useful in discriminating between cancerous and non-cancerous breast tissue. 2,6 This parameter is defined as γ ¼ ð1 − g 1 Þ∕ð1 − g 2 Þ, where g 1 and g 2 are the first and second Legendre moments of the phase function, respectively. (g 1 is also known as the anisotropy factor and is often written simply as g.) 23 The parameter γ characterizes the probability that a photon will undergo large-angle scattering during a scattering event.
Physiologically, research has shown that γ is related to the size-scale ratio of the particles in the tissue, which is the ratio of the tissue's largest particles to its smallest particles. 24 The diffuse photons measured by the aforementioned state-of-the-art SFDI methods are found in the 0 to 0.2 mm −1 spatial frequency range. These photons have traveled farther from their source of injection and deeper into the tissue (up to 24 mm deep in the absence of absorption) 25 before reflecting out of the tissue. As such, the individual angles they scatter with as they traverse through the tissue are averaged out over the large area of tissue through which they travel. Because of this averaging, reflectance in the diffuse domain is virtually insensitive to γ. 26 In order to capture information about γ, SFDI reflectance spectra must instead be measured and examined in the sub-diffuse domain, which typically comprises spatial frequencies of 0.5 mm −1 or greater for biological tissues. 26 SFDI that is inclusive of the sub-diffuse spatial frequency range is known as sub-diffuse spatial frequency domain imaging (sd-SFDI). Reflected photons found in the sub-diffuse range have traveled shallower into the tissue (<0.4 mm) 25 and undergone a smaller number of scattering events, making them more indicative of their initial scattering events. As such, they are more sensitive to γ and provide more information about the tissue at shallower depths. These photons are also sensitive to μ 0 s , but insensitive to μ a . 26 Recent work has explored capturing sd-SFDI images and analyzing the data to find γ and μ 0 s , 2,6,26 known as sub-diffuse optical properties. However, the great advancements that have been made to bring SFDI to real-time speeds in the diffuse domain have not been translated over to measuring and analyzing these spectra in the sub-diffuse domain. Improving the speed of sd-SFDI could open up a new opportunity in cancer treatment. Currently, there are two main methods by which surgeons excise cancerous tissue from a patient. The first is wide excision [ Fig. 1(a)], a procedure in which the surgeon removes, in addition to the tumor, a defined margin of normal tissue around the tumor to ensure that all cancerous tissue is captured. At the end of this procedure, the removed tissue undergoes histological analysis to confirm negative margins. 27 The other method is tissue-conserving surgery [ Fig. 1(b)], e.g., Mohs micrographic surgery. During this process, a small amount of tissue is removed, and the topmost surface of this resected tissue is analyzed histologically for negative margins. If cancer is present at the evaluated margins, the process is repeated until all of the cancerous tissue is removed and a negative margin is confirmed. Because of the repeated histological assessments required, this procedure is time and labor intensive. 28,29 The advent of wide-field sd-SFDI that identifies quantitative cancer biomarkers could allow for an improved method: image guided surgery [ Fig. 1(c)]. Image guidance could allow the surgeon to remove minimal amounts of tissue without the need for repeated histological analysis, combining the speed of wide excision with the mitigated normal tissue loss of tissue-conserving surgery.
In image guided surgery, resected tissue could be diagnosed with imaging tools in place of histology. Alternatively, the excision site itself can be examined if the analysis method does not require removal of the tissue. In both cases, the target area of examination would be the topmost layer of tissue. The sub-mm penetration depth of sd-SFDI makes this imaging modality an attractive option for this scenario. Additionally, sd-SFDI is high resolution (<1 mm), it has a wide field of view (>1 cm), and it is a slide-free, stain free method. 30 Image-guided surgery with sd-SFDI would require real-time measurements. More specifically, sd-SFDI images would have to be acquired in real-time and optical property heatmaps would have to be rendered from the sd-SFDI images in real time. 21 This paper focuses on the latter challenge. Previous experiments have demonstrated rendering tissue optical property heatmaps from sd-SFDI images using least-squares minimization on a semi-empirical equation. 2,6,26 The challenges to using this method are two-fold. Primarily, the method requires 15 min for a 300 × 300 pixel image. 26 Even for scenarios that do not require real-time imaging, such a long rendering time could serve as a significant bottleneck for advancing the knowledge of subdiffuse optical properties in biological tissue, since such studies in the diffuse domain typically entail large datasets and images much larger than 300 × 300 pixels. 20,21 Secondarily, the semiempirical equation was developed using data simulated with the modified-Henyey Greenstein phase function approximation, 26 whose domain 31 does not fully encompass the range of optical properties seen in biological tissue. 2,21,26,[32][33][34][35] Fig. 1 Flow diagram of existing and potential methods for removing cancerous tissue. All methods end in a final histological analysis to confirm all cancerous tissue has been removed. (a) In wide excision, a defined margin of normal tissue is removed along with the tumor in order to capture all cancerous tissue in one step. (b) In tissue-conserving surgery (e.g., Mohs micrographic surgery), small amounts of tissue are removed and then examined with histology repeatedly until all cancer is removed. (c) A potential method could use image guided surgery as a combined removal and examination step, resulting in a process that both captures all cancerous tissue in one step and mitigates normal tissue loss.
Recent work by Naglič et al. demonstrated a deep learning method for rendering sub-diffuse optical property heatmaps from sd-SFDI images that is capable of real-time speeds. 36,37 However, while this method was shown to be successful for simulated data, this method has never been calibrated for and applied to physical experimental data. Moreover, Naglič et al.'s model was trained on data simulated using the Gegenbauer kernel (GK) phase function approximation, whose domain 38 also does not cover the full range of optical properties seen in biological tissue. 2,21,26,[32][33][34][35] Another issue with sd-SFDI is that information on sub-diffuse optical properties is limited, and its use in cancer diagnostics has only been studied for breast tissue. 2,6 One avenue where sd-SFDI guided surgery could be particularly useful is in the field of skin cancer, where tissueconserving surgery is often performed using Mohs micrographic surgery. 28 Moreover, with extremely rare exceptions, skin cancer does not spread any other way but contiguously. 39 Because of this, if 100% of the excision surface is analyzed histopathologically, such as is done in Mohs surgery, 29 and no signs of cancer are found, no additional or deeper tissue analysis is needed to declare a negative margin. 28,29 This is in contrast to wide excision for skin cancer, which requires 4 to 6 mm of additional tissue to be removed and examined in order to declare a negative margin. 40 Not only does this advantage allow Mohs surgery to conserve more tissue, it also makes Mohs a prime candidate for use with sub-diffuse SFDI, as it makes Mohs amenable to the limited penetration depth of sd-SFDI. Previous research has shown that diffuse SFDI could be useful for discriminating between cancerous and normal skin tissue, [41][42][43][44] but little is known about the sub-diffuse optical properties of these tissue types.
This project is the first demonstration of real-time rendering of sub-diffuse optical property heatmaps from experimentally acquired sd-SFDI data. Moreover, we report preliminary data on the sub-diffuse optical properties of skin tissue subtypes. This work lays the foundation for implementing real-time applications of sd-SFDI, such as image-guided surgery.

Materials and Methods
In order to demonstrate real-time rendering of sub-diffuse optical property maps from sd-SFDI images, we acquired sd-SFDI images of two types of samples: tissue-simulating phantoms and human skin tissue samples. We imaged these samples across three wavelengths using a custom-built SFDI system. For each sample and wavelength, we collected 2050 × 2048 pixel demodulated reflectance images at 26 separate spatial frequencies, resulting in a [2050 × 2048 pixels × 26 spatial frequencies] image cube per sample per wavelength. We then developed a Monte Carlo model with novel phase function sampling and used it to generate a training and validation dataset of simulated [1 × 26 spatial frequencies] sd-SFDI reflectance spectra (RðfÞ). This can be thought of as simulating individual pixels of sd-SFDI image cubes and collecting them into datasets. We also simulated a reference phantom spectrum for calibration R ref ðfÞ for each wavelength. We used the training dataset to iteratively train an artificial neural network (ANN) to predict sub-diffuse optical properties from sd-SFDI measurements. The validation dataset was used to validate that the model's accuracy was improving during training iterations as to avoid overfitting the model to the training dataset.
We also fit a semi-empirical equation, developed by Kanick et al., 26 for benchmark comparisons. We tested the accuracy and speed of our ANN using the phantom measurements and compared the ANN's performance to that of Kanick's least-squared-minimization technique. 26 We then applied the model to cancerous skin tissue samples and visualized the differences in optical properties between cancerous skin tissue and other skin tissue subtypes.

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 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 mm −1 . The system would otherwise be limited to only ∼0.4 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 45 and made measurements less sensitive to variations in the samples' height. 46 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. 46 A 5-megapixel monochrome CCD camera (piA2400-17gm Basler, Ahrensburg, Germany) was used to collect images of the light reflected from the sample. 47 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. 11 This process was repeated for 26 spatial frequencies (0.01, 0.025 mm −1 , and 24 spatial frequencies evenly spaced between 0.05 and 1.25 mm −1 ) and 3 different wavelengths (450, 530, and 620 nm), resulting in 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.
an uncalibrated [2050 × 2048 pixels × 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 ref ðfÞ per pixel per wavelength. Each pixel of each uncalibrated image cube was then calibrated using the equation 11 resulting in a calibrated image cube for each wavelength.
The spatial frequencies used in this study include 6 frequencies in the diffuse domain (f ≤ 0.2 mm −1 ), 5 frequencies in the intermediate domain (0.2 mm −1 < f < 0.5 mm −1 ), and 15 spatial frequencies in the sub-diffuse domain (f ≥ 0.5 mm −1 ). The spatial frequencies in the diffuse domain were included as they improved estimation of μ 0 s (Supplement 1 in the Supplementary Material), which affects RðfÞ in both the sub-diffuse and diffuse domain. 2 This also meant that optical property estimation had to account for possible variations in μ a , which affects RðfÞ in the diffuse domain. 11 The upper bound of 1.25 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. 48 Optionally, additional demodulation processing 11 can be applied to the images to include the DC frequency, f ¼ 0 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 mm −1 which are already included. As such, we did not include images at f ¼ 0 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, 37 implemented with Keras and configured to work on OpenCL. 36 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 μ 0 s and γ, 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 μ 0 s and γ heatmaps at all three wavelengths.

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. 49 and Alerstam et al. 50 was modified to include the parameter γ. 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. 32 In order to allow for simulation across a wide range of the phase function parameters γ and g 1 , our look-up table is generated with the modified Henyey-Greenstein 23 (MHG) function when γ and g 1 fall within MHG's valid range (γ ≤ 1 þ g 1 31 ) and is otherwise generated with the GK. 51 We approximated GK's valid range to be γ > 1 þ 0.6g 1 and γ < 3 g 1 38 (Supplement 1 in the Supplementary Material), so our method cannot sample phase functions with γ ≥ 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][33][34][35] Our method uses a novel mapping which relates the γ and g 1 values to the GK parameters 38 (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 52 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. 32 This hybrid MHG-GK phase function sampling method allows for simulations of RðfÞ across a wide range of phase function parameters encompassing the γ and g 1 values seen in both our phantoms and in biological tissue. 2,21,26,[32][33][34][35] The method's ability to sample phase functions with a specified γ and g 1 combination, enabled through our GK mapping and MHG's analytical relation to γ and g 1 , 23 allows for the simulation of a reference phantom's spectrum R 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. 11 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 μ a values was included to account for variations in μ a values when estimating μ 0 s . 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 μ 0 s and γ. 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., γ ≥ 3 g 1 ) were discarded during simulation, shrinking the total number of samples in the training and validation dataset.
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, 35 and RðfÞ remains relatively invariant to changes in g 1 at g 1 values >0.8. 54 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 μ 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 μ 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, 26 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, 25 which range from around 1 to 2 mm in thickness. Fortunately, the impact this has on our model's μ 0 s 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 μ 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 μ a . Alternatively, a future study could use thicker tissue samples. Table 1 Parameters used for training, validation, and experimental datasets. Formats of "½x; yn", signify n elements in the range of ½x; y.

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 μm (SD: 0.01 μm), d 2 ¼ 0.19 μm (SD: 0.01 μm), and d 3 ¼ 0.99 μm (SD: 0.03 μm). We prepared two dilutions of each bead distribution stock solution in distilled water: one to achieve a value of μ 0 s ¼ 2 mm −1 at 530 nm and the other to achieve a value of μ 0 s ¼ 2 mm −1 at 530 nm, as calculated using Mie theory. 55 This resulted in six bead phantoms which covered γ values ranging from 0.97 to 2.17 and μ 0 s values ranging from 1.11 to 5.39 mm −1 to cover the ranges of values previously reported for various biological tissues. 2,21,26,[32][33][34] The optical properties of the phantoms at each wavelength were quantified using Mie Theory. 55 All phantoms had negligible absorption (assumed to be 0.001 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. 57 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 ref ðfÞ for each pixel and wavelength. When simulating this phantom to find R ref ðfÞ at each wavelength, μ a was set to 0.001 mm −1 to represent negligible absorption, as a value of 0 caused problems with the simulator.

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. 58 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. 56 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.

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 jεj ¼ 1∕n Ã P n i ðy i −ŷ i Þ∕y i . The model was then applied to the SFDI image cubes of the phantoms to render μ 0 s and γ heatmaps of the phantoms. The SFDI image cubes of the skin tissue were then processed using the trained ANN to render heatmaps of μ 0 s and γ. The skin tissue was assumed to have an anisotropy factor of g 1 ¼ 0.9. 45 Software timing tools were used to measure the frame rate at which the model was able to render these heatmaps.
In tandem, μ 0 s and γ predictions were rendered using the semi-empirical non-linear fitting method first reported by Kanick et al. 26 3 Results

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, 57 but it is often used to filter out specular reflections. 45 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). 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 μ 0 s and 3.6% for γ.
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 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 μ 0 s predictions and γ 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 jεj ¼ 0.77% for μ 0 s and jεj ¼ 0.94% for γ, on par with the accuracy that Naglič et al.'s model achieved over a narrower range of γ values. 36 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 × 2048 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. 14 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 nonlinear fit method processed a single pixel with an average speed 0.013 s∕pixel, while the ANN accomplished this in only 0.00187 s∕pixel. Thus even without parallelization, the ANN still outperforms the non-linear fit method by an order of magnitude.

Tissue Type Comparison
Some highlighted results of the tissue experiments can be seen in Fig. 4. We focus on μ 0 s and γ 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, μ 0 s heatmaps, and γ heatmaps for each tissue type. Areas of contrast can be seen in the heatmaps that align with the boundaries of the tissue subtypes.
Adipose tissue exhibits high μ 0 s , likely due to connective tissue septa which are particularly dense in subcutaneous fat in areas of the face and head, 41 from which 3 of our 4 samples were excised. Adipose tissue also exhibits high γ, indicating that the large adipocytes result in a high ratio of the adipose tissue's largest particles to smallest particles. 24 Dermis exhibits a relatively high μ 0 s and a low γ, which may be caused by its collagen fibers and relatively small constituent particles. 59 BCC tends to have a lower μ 0 s than normal skin tissue subtypes, echoing trends seen in the previous research. 41,43,44 It also exhibits low γ values, which may be associated with its high density of small nuclei. 2 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 μ 0 s values of BCC show a fair amount of overlap with those of adipose and dermis, the inner quartiles of μ 0 s values for BCC show almost no overlap with those of adipose and dermis. Meanwhile, γ 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 γ is a more reliable feature than μ 0 s 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 μ 0 s higher in adipose tissue than in BCC. However, this study also found that adipose tissue in skin excised from the back exhibited low μ 0 s values, on par with those of BCC, due to its thin connective tissue. 41 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 μ 0 s values. However, this tissue still exhibited high γ values similar to what is seen in our results. 2 These results indicate that γ could be more consistently higher in BCC than in adipose tissue, whereas μ 0 s 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 μ 0 s value versus the mean γ 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 mm −1 ) and a sub-diffuse frequency (f ¼ 0.7 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 μ 0 s versus γ 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. 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.

Discussions and Conclusions
These results represent the first time that sub-diffuse optical property heatmaps have been rendered from experimental sd-SFDI images in real time. This advancement, if combined with realtime sd-SFDI acquisition, could bring wide-field measurement of sub-diffuse optical properties to real-time speeds relevant for clinical applications (≥1 fps). 14 Real-time acquisition could be enabled by using state-of-the-art acquisition and demodulation techniques, 14,17 which can be readily used for sd-SFDI, and using fewer spatial frequencies in the ANN. Work by Naglič et al. demonstrates that sub-diffuse optical property heatmaps can be accurately rendered with ANNs using as few as five spatial frequencies: (0.0, 0.1, 0.2, 0.4, and 0.8 mm −1 ). 36 This heatmap rendering method can also be used for other applications outside of cancer detection, such as measuring burn wound severity 3 or agriculture applications such as detecting subsurface bruising in apples. 22,60 While the range of properties for which the model can produce accurate predictions is limited by the range of properties it is trained on, the model can easily be retrained on new datasets generated with the sd-SFDI Monte Carlo algorithm for any possible applications requiring a wider range of optical properties.
Our unique GK parameter mapping and combined MHG-GK phase function sampling method allow for targeted simulation of sub-diffuse reflectance spectra over a wider range of γ and g 1 values than previously possible. This sampling method enabled the application of the ANN to physical experimental sd-SFDI measurements. However, the method could also be used to simulate sub-diffuse spectra for other purposes, such as improving sub-diffuse optical property estimation from spatially resolved reflectance measurements. 61 Examination of the sub-diffuse optical properties of the tissue samples revealed intriguing differences between the properties of the cancerous skin tissue regions and those of the normal tissue regions explored in this study. Moreover, the observed trends are congruent with the limited existing information on sub-diffuse optical properties of biological tissue in general as well as on skin tissue specifically. 2,26,[41][42][43][44] These preliminary data are exciting, but measurements of a large number of tissue samples from several different subjects are needed to identify differences in the sub-diffuse optical properties of skin tissue subtypes that are both statistically and clinically significant.
Future studies can also expand this work to other skin tissue subtypes, such as epidermis and hair follicles. 41 Moreover, this work can be expanded to additional tissue types, with some caveats. Skin cancer is an ideal candidate for this work. There is already a well-established and widely used tissue-conserving surgery method, Mohs, in which examination of surface margins is standard, and the existence of remaining cancer can be confirmed or denied by examining the topmost layer alone. 28,29 However, the limited penetration depth of sd-SFDI is less applicable to cancers for which examination of the surface margin alone is insufficient to confirm negative margins, such as breast cancer. Breast cancer surgery requires examination of deeper tissue (i.e., up to 2 mm below the tissue surface for ductal carcinoma in situ) to obtain clinically accepted margin clearances. 62 Additionally, this work may be more difficult to incorporate for cancers for which tissue-conserving surgery is generally not considered a gold standard, such as with testicular cancer. 63 There are several other potential future studies which could be very useful. One would be to run a similar experiment using different wavelengths of light. This study used wavelengths in the visible range. These wavelengths penetrate tissue less deeply, 43 making them useful for taking the shallow measurements needed to examine surface margins during image guided surgery. However, skin tissue subtypes have been shown to have better separation of μ 0 s at longer wavelengths, 41 and there could exist a similar trend for γ. This possible trade-off would be worth exploring. Another future study which could improve these results further would be to incorporate profilometry correction, a technique which has already seen wide use in SFDI research. 21,46 All of our sd-SFDI measurements were captured with the angle of projection oriented orthogonally to the direction of the sine waves. Taking additional measurements where the angle of incidence is instead oriented parallel to the direction of the sine waves provides measurements that are highly sensitive to changes in the height of the sample and can be used to automatically correct measurements to account for these changes. 46 An additional reason to explore using sd-SFDI with multiple pattern orientations is images captured this way can be used to determine collagen fiber orientation, 64 information which is difficult to determine at wide-scale with other imaging techniques. 65 The structure of collagen fibers in skin tissue directly relates to the tissue's subtype, 41 and combining this information with the tissue's sub-diffuse optical properties could result in a powerful diagnostic feature set.
A possible limitation that we did not study is whether or not the freezing and thawing of the samples has any impact on the results of the optical property measurements. All samples were processed the same way, so any consistent effects from this would be consistent across all samples. However, it is possible that imaging the samples before freezing would yield different values of the optical properties. This could be problematic if that difference is significant, as the proposed use case of image-guided surgery does not involve freezing and thawing the tissue. There is also the possibility that optical properties found when imaging tissue that has not been frozen and thawed could be more accurate, or, more importantly, more useful in discriminating cancerous tissue. For these reasons, it would be useful to conduct a study on the impact of freezing and thawing a sample before determining its optical properties.
This study demonstrates a fundamental step forward in measuring sub-diffuse optical properties across a wide-field in real time and has already led to further projects in real-time sd-SFDI. 66 It also motivates further exploration into the potential of these properties in differentiating between cancerous and normal tissue subtypes. Altogether, it provides a foundation for the prospective goal of ultimately incorporating sd-SFDI into real-time medical applications, such as image-guided surgery.

Disclosures
The authors declare no conflicts of interest.
research interests include supervised machine learning and deep learning, and his graduate research applied these interests to sub-diffuse spatial frequency domain imaging.
Xu Feng received her PhD in biomedical engineering from The University of Texas at Austin and her MS degree in optical engineering from Tsinghua University. She is currently a postdoctoral fellow at Harvard Medical School and Massachusetts General Hospital. Her research interests include the development and translation of biomedical optical imaging techniques.
Yao Zhang received her BS degree in biomedical engineering from Shanghai Jiao Tong University in 2015 and her PhD in biomedical engineering from The University of Texas at Austin in Spring 2020 and works at Expedia Group as a data scientist. She is interested in translational research and using computational models, specifically machine learning models to solve problems.
Mia K. Markey is a professor of biomedical engineering at The University of Texas at Austin and an adjunct professor of imaging physics at the University of Texas MD Anderson Cancer Center. She is a fellow of both the American Association for the Advancement of Science and American Institute for Medical and Biological Engineering and is a senior member of both the IEEE and SPIE.
Biographies of the other authors are not available.