Optical sampling depth in the spatial frequency domain

Abstract. We present a Monte Carlo (MC) method to determine depth-dependent probability distributions of photon visitation and detection for optical reflectance measurements performed in the spatial frequency domain (SFD). These distributions are formed using an MC simulation for radiative transport that utilizes a photon packet weighting procedure consistent with the two-dimensional spatial Fourier transform of the radiative transport equation. This method enables the development of quantitative metrics for SFD optical sampling depth in layered tissue and its dependence on both tissue optical properties and spatial frequency. We validate the computed depth-dependent probability distributions using SFD measurements in a layered phantom system with a highly scattering top layer of variable thickness supported by a highly absorbing base layer. We utilize our method to establish the spatial frequency-dependent optical sampling depth for a number of tissue types and also provide a general tool to determine such depths for tissues of arbitrary optical properties.


Introduction
The use of spatial frequency domain (SFD) methods for diffuse optical imaging of biological tissues has gained significant traction in the biophotonics community since its introduction in 1998. 1 SFD methods combine a measurement of spatially modulated reflectance at multiple spatial frequencies with light transport models to determine optical and physiological properties of the tissue in question. The power of SFD methods is demonstrated most notably in spatial frequency domain imaging (SFDI), where such measurements are made for every pixel in a wide-field image. Such images provide functional mappings of optical and physiological properties with submillimeter detail 2 and temporal resolution limited to the spatial pattern projection times. 3 Diffuse optical methods utilize measurements made under multiple illumination/detection configurations whether it be multiple source-detector separations, delay times, and/or spatial/ temporal modulation frequencies. These configurations inherently collect photons that have penetrated different tissue volumes. Understanding the spatial regions that detected photons have sampled, and their sensitivity to source-detector configuration, is crucial in many situations, including the assignment of optical properties to a given tissue volume, accounting for the effects of tissue heterogeneities, measuring layered tissues, and the performance of image reconstruction. While this problem has been extensively studied for spatially/temporally resolved and temporal frequency domain methods, [4][5][6][7][8][9][10][11][12][13][14][15][16] extensive quantitative assessments of the optical sampling depths relevant to SFD methods are not prevalent in the literature.
Knowledge of tissue depths sampled by SFDI is important for understanding and contextualizing measurements of layered and heterogeneous tissue, and is relevant to all clinical and preclinical applications of SFDI described in the literature to date. For example, SFDI has been investigated for multiple applications in human skin, including burns, 17 reconstructive skin flaps, 18 and skin malignancies. 19 The layered structure of skin, which includes the largely avascular superficial epidermis and deeper vascularized papillary and reticular dermis and subcutaneous adipose tissue, requires knowledge of which layers are being probed in order to avoid misleading or irrelevant measurements. For example, in applications involving skin burns, the thickness of affected tissue dictates the treatment protocol, highlighting the importance of understanding the depth of tissue probed. 20 SFDI is also being explored for deeper tissue applications, including the measurement of human breast tumors, 21 where it is essential to understand the penetration of collected photons in order to evaluate the maximum depth of measurable tumor contrast. There has been a growing interest in subdiffusive SFDI, which makes use of high spatial frequencies (>0.2 mm −1 ) and has been used in applications such as tumor margin detection. 22 Again, quantitation of SFD sampling depth at these higher frequencies will allow users of this technique to determine the probed tissue thickness of resected specimens, which is relevant for interpretation of such results relative to the specific guidelines for clear margins for different tumor types. 23 In the preclinical setting, SFDI has been used for small animal tumor imaging to better understand cancer treatment response and resistance. 24 Photon sampling depth is essential for understanding what portion of the collected signal is due to superficial skin versus subcutaneous tumor. Finally, multimodality imaging, in which SFDI is combined with optical sectioning techniques such as multiphoton microscopy or optical coherence tomography, would benefit from knowledge of SFD photon sampling depth to ensure data from each modality is sampling similar tissue volumes. 25 Prior efforts to analyze SFD sampling depth utilized an approximate solution to the radiative transport equation (RTE) based on a spherical harmonic expansion 26 with errors dependent upon the order of the expansion. Here, we formulate a method to analyze SFD sampling depth by developing a Monte Carlo (MC) random variable that can be rigorously derived from the RTE. In this method, the error associated with the SFD sampling depth estimates is solely dependent on the number of photon packets launched. We provide both experimental validation for these results and sampling depth statistics for a range of optical properties in a manner that allows for rapid and accurate estimation for any candidate tissue.
We describe four steps toward the use of MC simulations to estimate SFD sampling depths. First, we provide a method to determine optical sampling depth statistics for SFD measurements using MC simulations. Second, we validate our simulated results by performing a series of experimental measurements on fabricated phantoms. Third, we apply our method to determine sampling depths for real tissue types. Finally, we provide tabular data and MATLAB code to generate sampling depth results for any tissue type. The archived version of the code can be freely accessed and executed through Code Ocean: https://codeocean .com/capsule/97500821-e8c9-456c-9a41-3926677bfc79/code.

Methods
We analyze the sampling depth of SFD measurements using MC simulations that determine the sampling probability of detected photon packets for specific depths. For this analysis, we utilize the modified shortcut method (MSM) 27 that directly performs an MC simulation of the two-dimensional (2-D) spatial Fourier transform of the RTE. Thus, the MC simulations are performed directly in the SFD and not subject to inaccuracies that can result from computation of the discrete Fourier transform of spatially resolved reflectance simulation data. Our computational model utilizes a representation that segments the tissue into z-plane surfaces at uniform depth increments and determines the subset of detected photon packets that traverse each depth. While the computational system that we consider is spatially homogeneous, the approach that we introduce is general and applicable to any layered tissue geometry.
To validate the computational predictions, we utilize a twolayer phantom system with a highly scattering layer of water, nigrosin, and TiO 2 particles placed over a highly absorbing layer of nigrosin-doped agar gel. We provide results of SFD measurements performed on this two-layer phantom system with differing thicknesses of the top layer for comparison with the MC simulation results.
Below, we establish a rigorous metric for depth-dependent probability of photon visitation and detection P V∩D ðzÞ and the use of an MC simulation for its computation. This depthdependent probability distribution is then related directly to the measured SFD reflectance and used to define various metrics for SFD optical sampling depth. We then discuss the details of our SFD validation measurements.

Monte Carlo Probability of Visitation and Detection P V∩D
We wish to characterize the spatial distribution of only those photons that are launched by a specified source and subsequently captured by a detector of interest. We accomplish this by computing a "contributon" response function that was first developed in the nuclear engineering community for the solution of deep-penetration nuclear transmission problems using MC methods. [28][29][30][31][32][33][34] In the transmission problems, a surface between the source and detector is first specified. A forward simulation for photon transport from the source is then matched with an adjoint simulation of photon transport from the detector over the midway surface using the contributon response function. In a reflectance geometry, we use this idea to determine the probability that detected photon packets have visited a depth d within the tissue by defining the surface at the selected depth as our midway surface. We determine the probability that photon packets from the source "visit" d, PðVÞ, and then determine the probability that they will subsequently be detected, PðDjVÞ. Bayes theorem 35 is used to determine the probability of visitation and detection, P V∩D E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 5 1 4 P V∩D ¼ PðVÞPðDjVÞ: (1) One method to determine only those photon packets that originate from the source, travel to the midway surface, and subsequently arrive at the detector, would be to match the radiance determined by a forward simulation from the source with an adjoint simulation from the detector over the midway surface. This approach works well when the source and detector are "small" relative to the midway surface. 14,36 In the biomedical optics community, the use of such coupled forwardadjoint simulations has been used to address fluorescence excitation and detection, 37 "photon hitting density" maps, 6 and tomographic sensitivity analysis for spatially resolved reflectance. 14,36,38 In this work, we produce P V∩D distributions for SFD measurements using a single conventional MC simulation. We take this approach because MC simulations of SFD methods utilize a "small" source and "large" detector, i.e., light is injected into the medium at a single point while detection occurs at all locations on the tissue surface. In such a scenario, a conventional MC approach provides better computational efficiency relative to coupled forward-adjoint methods. 38 Note that the SFD situation is unique relative to other diffuse optical methods where both "small" sources and "small" detectors are typically used. Specifically, we create depth-dependent SFD P V∩D ðzÞ distributions by defining parallel x − y planes, placed at regular intervals of 0.01 mm below the tissue surface, as our midway surfaces. We choose this interval size to provide submillimeter resolution for our probability distributions. We then perform an MC simulation that tracks those photon packets that both visit these depths and are subsequently detected. Figure 1 shows a schematic of the tissue segmented into layers using planes located at depths z 0 ¼ 0; z 1 ; z 2 ; · · · separated by height Δz. These planes define midway surfaces at depths z i at which P V∩D ðz i Þ is determined. We calculate P V∩D ðzÞ by tracking each photon packet from the source to these various midway surfaces and from those surfaces to the detector. In doing so, these calculations inherently combine the probability of a photon packet visiting the midway surface after being launched by the source with the probability of detection after visiting the midway surface. The P V∩D ðz i Þ tallies are computed using an MC simulation of radiative transport by determining the midway surfaces visited by each detected photon packet. As mentioned earlier, we utilize an MC simulation using the MSM 27 that tallies the photon packet weights directly in the SFD. A narrow, collimated beam normally incident on the tissue surface is used as the source within the simulation and the photon packets are transported using conventional MC propagation. 39,40 The complex weights specified by the MSM for the specified spatial frequency provide both the amplitude change and phase shifts that result in photon packet propagation from source to detector. When utilizing illumination modulated along the x-axis at spatial frequency f x , the SFD complex tally is 27 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 4 4 4 where ξ is a random variable that represents the tally or detected weight for each photon packet, W is the weight of an individual photon packet as determined by discrete absorption weighting, 41 x 0 is the location along the x-axis where the photon packet enters the tissue sample, and x f is the exiting x-axis location immediately prior to detection. Note that Eq. (2) can be used in layered tissue systems. The optical property changes in each layer will be implicitly captured by W and the location of x f . For each detected photon packet, the SFD complex random variable [Eq. (2)] is tallied for each midway surface that the photon packet crossed and notated as ξðz i Þ. For example, if the photon packet reached a maximum depth z max prior to detection, the photon packet crossed all midway surfaces residing at depths less than z max . For the photon packet trajectory shown in Fig. 1, the detected photon packet weight is tallied to the midway surface z 3 as well as all midway surfaces residing at shallower depths because the photon packet also crossed those surfaces. Once N photon packet trajectories are simulated, the expected value of ξ for depth z i , E½ξðz i Þ produces the probability that the trajectory of a detected photon packet crossed depth z i , P V∩D ðz ¼ z i Þ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; where N is the number of photon packets launched. By determining ξðz i Þ for all the z i depths under consideration, we form the depth-dependent probability distribution of photon packet visitation and detection ¼ P V∩D ðzÞ.

Monte Carlo Maximum Depth of Penetration P z max
The P V∩D ðzÞ distribution over all depths z does not result in a directly measurable quantity because photon packets that contribute to the P V∩D tally at a surface z i also contribute to P V∩D tallies at all locations shallower than z i . However, we can use P V∩D ðzÞ to derive a distribution that isolates the contribution of each detected photon packet to a single bin corresponding to the maximum depth visited by that photon packet trajectory.
In doing so, we ensure that each detected photon packet is tallied to only a single depth bin within the tissue. This distribution is formed by taking differences of the P V∩D ðzÞ tally in successive bins E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 6 0 2 We call P z max ðzÞ the "z max " distribution which isolates those detected photon packets that crossed into depth z i but did not cross into z iþ1 . Because each photon packet is tallied only once at its maximum depth of propagation "z max ," P z max ðzÞ has the property that its integral over all depths z recovers all the weight of all the detected photon packets, which is equivalent to the total diffuse reflectance, R d E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 4 9 4 If the upper limit of integration on the right-hand side of the above equation is taken instead to some finite depth d, the result would tally only those photon packets that contribute to the reflectance and whose trajectories were restricted to depths less than or equal to d. We notate this as P z max ðz ≤ dÞ or E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 3 9 5 This construct will be used to validate our MC simulation results with our SFD measurements. Division of P z max ðz ≤ dÞ by R d produces a probability distribution function that describes the fraction ðXÞ of the detected light that visited tissue depths d or less E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 2 9 7 By setting the value of X in the above equation to a given value, say 0.5 (or 50%), we can calculate the maximum tissue depth d 50 from which 50% of the detected reflectance emanates.

Experimental Validation
We used experimental SFD measurements to validate our computational results. These measurements were taken in a twolayer phantom designed to determine the optical sampling depth as a function of spatial frequency. The two-layer phantom was housed in a container with (L × W × H) dimensions of 7.2 cm × 10.8 cm × 6.1 cm. The top layer of the phantom was a liquid composed of water, nigrosin, and TiO 2 particles with optical properties μ 0 s ∕μ a ¼ 100 and l Ã ¼ 1∕ðμ a þ μ 0 s Þ ¼ 2 mm at λ ¼ 731 nm. The bottom layer was a highly absorbing solid phantom composed of agar, water, and nigrosin, and occupied a total volume of 350 mL in the container. The top layer Journal of Biomedical Optics 071603-3 July 2019 • Vol. 24 (7) thickness d was varied from ½0À7.5l Ã . Figure 2 provides a schematic of the setup. We used the OxImager RS SFDI system (Modulated Imaging Inc., Irvine, California) to measure the two-layer phantom. The SFDI system utilizes crossed linear polarizers in front of the projection and detection lenses to minimize the effect of specular reflection and select for diffuse reflection. The projection field of view (FOV) was 20 cm × 15 cm and directed to the tissue surface at an angle of 15 deg relative to the surface normal. Detection was performed perpendicular to the surface of the phantom with an FOV of 8 cm × 6 cm (effective NA = 0.253). A vertical translation stage was used to support and adjust the height of the container with the two-layer phantom. After taking the first measurement of the highly absorbing solid phantom (i.e., d ¼ 0 mm), incremental amounts of the liquid phantom were added successively with a predetermined volume such that the top liquid layer thickness d above the solid phantom increased by 0.5 mm between each measurement following the d ¼ 0 mm measurement. The micrometer on the translation stage was used to lower the two-layer phantom system by 0.5 mm following each measurement such that the top surface of the two-layer phantom system remained at a constant image plane for all top layer thicknesses measured. This was done to avoid the need for height correction during data processing. All measurements were performed at a wavelength λ ¼ 731 nm with spatial frequencies f x ¼ 0, 0.0125, 0.025, 0.0375, 0.05, 0.0625, 0.075, 0.0875, 0.1, 0.125, 0.15, 0.175, 0.2, 0.25, and 0.3 mm −1 . At each spatial frequency (in this case, along one spatial dimension x), raw reflectance images at three different phases (0, 2π∕3, and 4π∕3 radians) were sequentially projected onto the phantom using a digital micromirror device. The resulting reflected light was imaged with a camera. The images were then demodulated to extract the amplitude envelope for each spatial frequency measurement using an established amplitude demodulation algorithm. 42,43 A separate reference measurement at the same spatial frequencies was made on a calibration phantom with known optical properties for calibration of the SFDI source intensity and instrument response. This calibration enables the measured reflectance to be converted to absolute reflectance. This is achieved by comparing the measured reflectance from the calibration phantom with its predicted diffuse reflectance from an MC-based forward model using the phantom's known optical properties. The region of interest (ROI) chosen for data analysis was centered in the detection FOV to avoid edge effects and measured 7.5 cm × 2 cm. The reported experimental data are average values taken over the ROI.

Results and Discussion
We first present P V∩D ðzÞ results from which we compute P z max ðzÞ. The P z max ðzÞ results will form the basis for (a) validating our computational method with experimental measurements and (b) determining metrics for optical sampling depth. We then present optical sampling depth results for a variety of tissue types based on literature reported optical properties.

Probability of Visitation and Detection P V∩D
We first consider a highly scattering tissue system with refractive index n ¼ 1.4 with optical properties providing ðμ 0 s ∕μ a Þ ¼ 100 and l Ã ¼ 1 mm and single-scattering anisotropy g ¼ 0.8. This corresponds to optical absorption and scattering coefficients μ a ¼ 0.00990099 mm −1 and μ s ¼ 4.95049505 mm −1 , respectively. In the MC simulation, we utilize a narrow collimated beam normally incident into the tissue. Using the method described in Sec. 2.1, we launched N ¼ 10 8 photon packets to obtain P V∩D ðzÞ with z bins incremented at 0.01 mm and a set of spatial frequencies spanning 0 and 0.5 mm −1 .  to the total diffuse reflectance for each spatial frequency, R d ðf x Þ. This is because every diffusely reflected photon packet will traverse this surface. The plots are shown with 1 − σ error bars which indicate that 68% of independent simulations will produce results, which lie within this interval. The f x ¼ 0 mm −1 plot has a value of 0.62 at z ¼ 0 mm representing total diffuse reflectance and decays to 4.7 × 10 −4 by depth z ¼ 20 mm, representing approximately a 3 order of magnitude reduction. The f x ¼ 0.5 mm −1 plot has a value of 0.03 at z ¼ 0 and decays to 3.6 × 10 −7 at z ¼ 20 mm, nearly a 6 order of magnitude reduction. These data illustrate the low-pass optical transport characteristics of scattering tissues. Figure 3(b) shows the P z max ðzÞ derived from P V∩D ðzÞ using Eq. (4). As described in Sec. 2.2, the benefit of transforming the P V∩D ðzÞ results to the max z i formulation P z max ðz ¼ z i Þ is that (a) P z max ðz ¼ z i Þ, when integrated over all z i , produces total diffuse reflectance R d and (b) it effectively tracks the maximum depth sampled by each photon packet and provides group statistics on sampling depth for each spatial frequency.

Experimental Validation
The experimental setup described in Sec. 2.3 consists of measurements taken from a two-layer phantom with a highly absorbing bottom layer placed at various depths d that extinguishes any photons that propagate to that depth. The resulting measured reflectance is composed of only photons that never reach depths z > d, i.e., the photons detected possess trajectories with a maximum z ≤ d. The analogous computational result is given by Eq. (6). Figure 4(a) shows the calibrated experimental measurements of diffuse reflectance versus top layer thickness d and plots of P z max ðz ≤ dÞ and (b) their difference. The plot shows a subset of the measured f x values for clarity. The depths d and the spatial frequencies have been normalized to l Ã ¼ 1 mm. For normalized spatial frequency f x l Ã ¼ 0, the plot rises from 0 and reaches a value of 0.595 for top layer thickness of 7.5 d∕l Ã . The plot rises monotonically as the top layer thickness increases because increasing numbers of photons fail to be extinguished by the bottom layer and are able to return to the surface to contribute to reflectance. As the spatial frequency increases, the measured diffuse reflectance flattens at a certain depth indicating that spatially modulated light for this frequency does not interrogate the tissue below that depth. For example, for f x l Ã ¼ 0.3, the spatially modulated reflectance rises from 0 and rises to 0.05 for top layer thickness of 1.0 d∕l Ã without further increases for larger top layer thicknesses. The absolute difference between the experimental measurements and the computational predictions ranges between [−0.012; 0.025].

Metrics for Optical Penetration Depth
We can determine the depth-dependent variation of the detected photon packets by determining the tissue depth that contains the complete photon packet trajectories corresponding to a certain fraction of the total measured diffuse reflectance at a given spatial frequency. For example, the spatial frequency dependence of sampling depth that encloses all the photon packet trajectories corresponding to only 50% of the detected reflectance can be determined by finding the value d 50 that results in a value of X ¼ 0.5 using Eq. (7). Similarly, by substituting alternate values for X ¼ 0.1, 0.25, 0.75, and 0.9 into Eq. (7), we computed depths d 10 , d 25 , d 75 , and d 90 that correspond to the tissue depths that enclose the photon packet trajectories responsible for 10%, 25%, 75%, and 90% of the detected reflectance, respectively. These are shown in Fig. 5 with specific numerical values provided in Table 3. The span of the gray rectangles [25 to 75]% and vertical-capped lines [10 to 90]% provides range of depths sampled by these portions of the detected reflectance for these optical properties. In Appendix B, we provide similar plots and tables for a several μ 0 s ∕μ a values to show how these spans vary with optical properties. These depth metrics are determined from the z max distribution and provide a direct correspondence between the maximum tissue depths sampled by portions of the detected reflectance. Unlike the photon hitting density, 6 these depths do not tally the pathlengths that the detected photon packet traverses within   the tissue volumes considered and do not map directly to absorption sensitivity. Rather they provide tissue depths that correspond to portions of the detected reflectance. These allow the SFD user to determine what fraction of the measured reflectance is available to interact with the tissue beyond a certain depth. For example, if a given tissue has a d 90 ¼ 5 mm, the user is certain that 90% of the measured reflectance is restricted to tissue depth ≤5 mm whereas only 10% of the measured reflectance has the opportunity to sample tissue depth >5 mm.

Sampling Depth Estimates for Real Tissue Types
The results provided in Fig. 5 can be generated for any set of optical properties. We apply our methodology to various tissue types to determine the depths that correspond to detection of 50%, 75%, and 90% of the total reflectance at a given spatial frequency. Table 1 shows the tissue types and wavelengths considered and the corresponding optical properties. [44][45][46][47] We performed independent MC simulations for each tissue type and wavelength pair to determine depth estimates for each. Figure 6 shows the predicted sampling depths versus spatial frequency for the various pairs of tissue type and wavelength considered. The tissue type, in which spatially modulated illumination can probe most deeply at the smaller spatial frequencies, is human breast at λ ¼ 851 and 731 nm. This is due to the higher μ 0 s ∕μ a properties for this tissue, 256 and 145, respectively, along with moderate transport mean-free path values of l Ã ¼ 0.89 mm and 1.04 mm, respectively. The median sampling depth of the human brain is 20% smaller than for human breast. While human brain tissue at λ ¼ 731 nm has the equivalent μ 0 s ∕μ a value as human breast at λ ¼ 851 nm, the transport meanfree path in human brain is only l Ã ¼ 0.76 mm as compared to l Ã ¼ 1.04 mm in human breast. While mouse skin has the lowest μ 0 s ∕μ a values of the four tissue types considered, which would suggest more superficial optical sampling depths, the transport mean-free paths are the largest of the tissues considered resulting in optical sampling depths that are only slightly more superficial than human brain tissue at low spatial frequencies. The μ 0 s ∕μ a values for human skin are not much larger than mouse skin but with much higher scattering properties resulting in the smallest l Ã values of the group and the most superficial optical sampling depths.
The spatial frequency dependence characteristics of these optical sampling depths are also of interest. For spatial frequencies larger than f x ¼ 0.1 mm −1 , differences in the optical sampling in human breast, human brain, and mouse skin are barely distinguishable. For human brain, the median depth values for

Determination of Optical Sampling Depths for Other Tissue Types
To enable the determination of optical sampling depths for any tissue type, we used our methodology to determine the variation of optical sampling depth with spatial frequency for a range of μ 0 s ∕μ a values between 1 and 1000 while keeping l Ã ¼ 1 mm fixed and assuming g ¼ 0.8 and n ¼ 1.4. The specific optical properties considered are listed in Table 2. Figure 7 shows the median depth of optical sampling as a function of μ 0 s ∕μ a . In this figure, both the median sampling depth and the spatial frequency of illumination are scaled relative to l Ã . The plots show that as the μ 0 s ∕μ a increases, so does the median depth. Moreover, at larger spatial frequencies, there is less sensitivity of the sampling depth to variations in μ 0 s ∕μ a . Optical sampling depths determined using the general optical properties can be composed into a lookup table, then scaled, and interpolated to provide depth statistics for an arbitrary tissue. Only knowledge of the tissue absorption and reduced scattering properties is needed. The details of how this is performed are given in Appendix B. The advantage of this lookup table method is that it dispenses with the need to execute an MC simulation for the specific tissue optical properties in question and enables rapid estimation of SFD sampling depths. We performed the lookup table method using the μ a and μ 0 s values of the real tissue optical properties listed in Table 1. Figure 8 shows the relative differences between the median depth d

Conclusions and Future Work
We have presented a transport-rigorous MC method to determine optical sampling depth statistics in the SFD. This method provides depth-dependent probability distributions of photon visitation and detection [P V∩D ðzÞ] for each spatial frequency within homogeneous or layered tissue. Our sampling depth predictions were validated experimentally using SFD measurements taken on a custom fabricated two-layer phantom system. Excellent agreement was obtained between these measurements and our MC predictions.
We applied our method to provide depth sampling statistics for a variety of tissue types at commonly used wavelengths. We nondimensionalized our results to create a 2-D lookup table to determine sampling depth statistics for any tissue given knowledge of the absorption and reduced scattering properties of the tissue. We provide this table and associated computer code to enable its use in the supplemental material.
Collectively, this work provides a rigorous methodology and convenient means to determine optical sampling depth in the SFD. Moving forward, we wish to analyze the effect of spatial-frequency-dependent variations in optical sampling depth on the extraction of optical properties when using SFD measurements at two or more spatial frequencies. 43 The use of SFD measurements at multiple spatial frequencies results in differential penetration depths leading to a partial volume effect for measurements taken in heterogeneous media. This effect is potentially reduced by choosing spatial frequencies proximal to each other, but this can compromise the ability to accurately extract optical properties. 48 Future work will aim to evaluate the tradeoffs between partial volume effects and optical property extraction errors when choosing spatial frequencies for specific applications. Table 2 General optical properties with l Ã ¼ 1 mm used for our sampling depth lookup table.  Appendix A: Real Tissue Data Details In Figs. 9-16, we provide detailed sampling depth metrics for the specific tissue types and wavelengths as listed in Table 1 with partial results shown in Fig. 6

Appendix B: General Tissue Data Supplemental Tables and Code
We generated a table of data for general tissue optical properties that can be used to estimate sampling depths for any arbitrary tissue i given knowledge of the μ a;i and μ 0 s;i of the candidate tissue. The data consist of 2-D lookup tables for the 10%, 25%, 50%, 75%, and 90% sampling depths (Table 3). Each table has μ 0 s ∕μ a along one axis and f x l Ã along the other. The μ 0 s ∕μ a values are those listed in To determine the sampling depth for a candidate tissue with optical properties μ a and μ 0 s , each table entries and axes are scaled appropriately and then a 2-D linear interpolation method, interp2 (MATLAB 2016b), is used to determine the sampling depth at the spatial frequency of interest. Specifically, to scale the table appropriately, the entries are multiplied by l Ã of the candidate tissue, and the f x l Ã axis is divided by l Ã of the candidate tissue. Then the μ 0 s ∕μ a of the candidate tissue and the spatial frequencies f x of interest are used to interpolate into each table to produce sampling depth estimates.
The tabular data and code to interpolate the data are given in the supplemental material.      (7)