Model for the diffuse reflectance in spatial frequency domain imaging

Abstract. Significance In spatial frequency domain imaging (SDFI), tissue is illuminated with sinusoidal intensity patterns at different spatial frequencies. For low spatial frequencies, the reflectance is diffuse and a model derived by Cuccia et al. (doi 10.1117/1.3088140) is commonly used to extract optical properties. An improved model resulting in more accurate optical property extraction could lead to improved diagnostic algorithms. Aim To develop a model that improves optical property extraction for the diffuse reflectance in SFDI compared to the model of Cuccia et al. Approach We derive two analytical models for the diffuse reflectance, starting from the theoretical radial reflectance R(ρ) for a pencil-beam illumination under the partial current boundary condition (PCBC) and the extended boundary condition (EBC). We compare both models and the model of Cuccia et al. to Monte Carlo simulations. Results The model based on the PCBC resulted in the lowest errors, improving median relative errors compared to the model of Cuccia et al. by 45% for the reflectance, 10% for the reduced scattering coefficient and 64% for the absorption coefficient. Conclusions For the diffuse reflectance in SFDI, the model based on the PCBC provides more accurate results than the currently used model by Cuccia et al.


Introduction
Spatial frequency domain imaging (SFDI, also known as modulated imaging or structured light imaging) is a technique that facilitates wide-field imaging of tissue optical properties 1,2 and has gained considerable interest in recent years for a range of applications, 3 such as tumor margin assessment in breast cancer 4 and assessment of diabetic feet. 5In SFDI, tissue is illuminated with sinusoidal intensity patterns and the reflected intensity is captured by a camera.Scattering and absorption of light by the tissue change the amplitude of the reflected intensity pattern (but not its spatial frequency).The amplitude of the reflected intensity patterns, M AC , depends on the tissue optical properties, the projected spatial frequency, and the properties of the optical system itself.By measuring the reflectance at two or more spatial frequencies, the tissue optical properties can be extracted from SFDI measurements.Compared to other wide-field optical techniques, such as conventional multi-and hyperspectral imaging, SFDI has the advantage that it can be performed using only a few wavelengths, improving acquisition time.Furthermore, by changing the projected spatial frequency, the interrogation depth of SFDI can be modified 6 to match the clinical application.
To extract optical properties from SFDI measurements, the reflected intensity is first demodulated.The amplitude modulation M AC for any tissue location x, and projected spatial frequency f x can be expressed as 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 ; 1 1 6 ; 6 9 9 tissue ðx; f x ; μ a ; μ s ; pðθÞÞ; (1) where I 0 is the projected intensity pattern, MTF SYS is the modulation transfer function of the system, and R tissue is the reflectance from the tissue which also depends on the tissue optical properties-the absorption coefficient μ a , the scattering coefficient μ s , and the phase function pðθÞ.The phase function describes the probability of scattering at a certain angle θ with respect to its previous direction.In general, two methods exist to demodulate the reflectance.For a detailed description, we refer to. 3 In the single-pixel demodulation method, three images are acquired with the same projected spatial frequency, but phase-shifted 2π∕3.From these three images, M AC is calculated through a simple analytical function.In the multipixel demodulation method, a single image is acquired and M AC is calculated by taking a Fourier transform of a line or the entire image.
To obtain R tissue from M AC , the modulation transfer function of the optical system has to be determined.MTF SYS is obtained through a calibration step where the M AC of a reference sample is measured for which the optical expected reflectance R ref is known, 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 ; 1 1 6 ; 5 1 3 To extract tissue optical properties from SFDI measurements of R tissue , two main approaches exist.The first approach is to use an analytical model based on physics that relates R tissue to the tissue optical properties. 2,7][10] In this manuscript, we will focus on the first approach, an analytical model based on physics.
Two analytical models exist for SFDI, the model of Cuccia et al. for the diffuse regime 2 and the model of Kanick et al. for the subdiffuse regime. 7The model of Kanick et al. is a semiempirical model based on MC simulations and the model of Cuccia et al. was derived from diffusion theory.Whether measurements are in the diffuse or subdiffuse regime depends on the projected spatial frequency and tissue optical properties.In this manuscript, we will focus on modeling the diffuse reflectance in SFDI.Diffusion theory is generally thought to be valid for f ≪ μ tr and μ 0 s ≫ μ a (where μ tr ¼ μ a þ μ 0 s ).Under the assumption of linearity of the medium (i.e., the frequency and phase of the modulated incident light are maintained in the modulated fluence rate in the tissue), Cuccia et al. obtained the following diffusion equation: 2 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 ; 1 1 6 ; 2 8 8 ∇ 2 z φðzÞ − μ 02 eff φðzÞ ¼ −3μ tr SðzÞ; where φðzÞ is the fluence rate at depth z, p is the effective attenuation coefficient that takes into account the influence of the projected spatial frequency , and D is the diffusion coefficient D ¼ 1∕ð3ðμ 0 s þ μ a ÞÞ, SðzÞ is a source term, and μ tr is the transport coefficient.The only difference with the diffusion equation for a uniform plane illumination (k ¼ 0) is in the scalar parameter μ 0 eff .Therefore, any approach to obtain the reflectance for a uniform plane illumination can be used for SFDI by replacing μ eff with μ 0 eff .For the diffuse reflectance for a uniform plane illumination, Cuccia et al. used the model of Svaasand et al. 11 for the diffuse reflectance for an infinitely wide illumination source and obtained the following model: 2 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 ; 1 1 6 ; 1 4 1 where a 0 ¼ μ 0 s ∕μ tr is the reduced albedo and A ¼ ð1 þ R eff Þ∕ð1 − R eff Þ describes the influence of the refractive index mismatch between the sample and the surrounding medium.Please note that a different expression for A was used by Cuccia et therefore Eq. ( 4) might seem different from their manuscript but the equation is only rewritten.In this paper, we propose a new model for the diffuse reflectance in SFDI and show that it reduces the median error in the estimated reflectance and extracted optical properties.

Modeling the Diffuse Reflectance in SFDI
As an alternative to directly solving the diffusion equation for a uniform plane illumination, a tissue "impulse response" can be obtained by computing the reflectance as a function of radial distance, RðρÞ, for illumination by an infinitely narrow pencil beam.The reflectance as a function of spatial frequency is then found by a 2D Fourier transform, which in the case of cylindrical symmetry is the same as computing the zeroth order Hankel transform: ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 5 9 7 where J 0 is the zeroth order Bessel function of the first kind.This approach is often used to translate MC simulations that generate RðρÞ to reflectance values R tissue ðkÞ for SFDI for a range of spatial frequencies. 2,7The same principle can equally well be applied to analytical solutions of RðρÞ.
In general, diffusion theory assumes infinite media.For semi-infinite media such as tissue, the refractive index mismatch between the sample and the medium above it results in a significant fraction of radiant energy being reflected back into the sample upon interaction with the boundary.Analytical expressions for RðρÞ can be obtained by imposing appropriate boundary conditions at the interface between the sample and the (non-scattering) medium above it.Analytical expressions for RðρÞ in the diffuse regime are available for different (yet equivalent) boundary conditions, i.e., the partial current boundary condition (PCBC) as proposed by Keijzer et al. 12 and the extended boundary condition (EBC) as proposed by Farrell et al. 13 For the PCBC, the irradiance at the boundary is set equal to the integral of the reflected radiance 12 and for the EBC, the fluence rate is set to zero at an extended boundary located at a distance z b outside the sample.There is no theoretical reason to prefer the PCBC over the EBC, the difference merely demonstrates a limitation of diffusion theory. 14The two boundary conditions result in the following expressions for the reflectance: 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 ; 1 1 6 ; 3 4 9 and 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 ; 1 1 6 ; 2 9 3 where A describes the influence of the mismatch between the refractive index of the sample, n i , and the surrounding medium, n e , which is equal to 15,16 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 2 0 8 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 1 4 6 where θ r is the angle of refraction.
The reflectance as a function of radial distance is then found as RðρÞ ¼ ∫ ∞ z 0 ¼0 Rðρ; z 0 ÞSðz 0 Þdz 0 , where S is the source term, for which we use a distributed source term SðzÞ ¼ a 0 μ tr expð−μ tr zÞ.To obtain our new models for the SFDI reflectance we integrate RðρÞ for each boundary condition over the radial coordinate ρ from 0 to ∞ and we replace μ eff with μ 0 eff .Performing the integrations for R PCBC ðρ; z 0 Þ and R EBC ðρ; z 0 Þ over z 0 and ρ we arrive at E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 6 6 0 and ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 5 9 5 3 Methods Since the model of Cuccia et al. was based on the PCBC, we first compare R Cuccia ðkÞ and R PCBC ðkÞ to the Hankel Transform of R PCBC ðρ; z 0 Þ.To calculate this Hankel Transform with Matlab we first integrated R PCBC ðρ; z 0 Þ over z 0 from 0 to 100 mm and r from 0 to 500 mm (increasing the integration limits did not change the results).Next, to compare the accuracy of the three different models (R Cuccia ðkÞ, R PCBC ðkÞ and R EBC ðkÞ) we performed MC simulations to obtain the reflectance versus radial distance, R MC ðρÞ.To calculate the reflectance measured by SFDI, R MC ðkÞ, we performed a Hankel Transform.For each model, we calculated the relative error (RE) in the reflectance with respect to R MC ðkÞ for all the simulations E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 4 1 0 We also determined the relative errors in the values of μ a and μ s 0 that would be obtained for each of the models if they would be used to fit measured reflectance values obtained with spatial frequencies of 0 and 0.5 mm −1 .

Monte Carlo simulations
We simulated a pencil beam illumination and collected photons versus radial distance from the source with a 0.001-mm bin size and 4 • 10 5 bins (regardless of their angle upon detection) and we performed a Hankel Transform to obtain the reflectance for a given spatial frequency.We simulated tissues with all combinations of μ a ¼ ½0.001; 0.005; 0.01; 0.05; 0.1 mm −1 , μ 0 s ¼ ½1; 5; 10; 20; 50 mm −1 , and two different phase functions.One set of simulations was done with a Henyey-Greenstein (HG) phase function with g 1 ¼ 0.9, and a second set of simulations was done with a two-term HG (TTHG) phase function since the majority of published phase function measurements are best described by a TTHG. 17We used a TTHG phase functionwith a scattering anisotropy g 1 of ∼0.83, using the following parameters: pðθÞ ¼ 0.45 • P HG ðg HG ¼ 0.95Þ þ 0.05 • p HG ðg HG ¼ −0.2Þ,where p HG denotes a regular HG phase function.For each set of optical properties, we used spatial frequencies such that μ 0 s f −1 ranged from 0.1 to 1000 with 20 equal steps on a log-scale for each value of μ s 0 .We simulated a refractive index of the tissue of 1.33 and the medium above it of 1.00 (A ¼ 2.515).For diffusion theory to hold, f ≪ μ tr 0 , so we excluded simulations where f > 1 5 • μ tr 0 .We performed each simulation three times and ensured we launched enough photons so that the standard deviation over the mean reflectance was <1%.To determine the accuracy of the models we compared these to the reflectance values averaged over these three simulations.
We first compared R Cuccia ðkÞ [Eq.( 4)] and R PCBC ðkÞ [Eq.(10)] to the Hankel Transform [Eq.( 5)] of R PCBC ðρ; z 0 Þ [Eq.( 6)], since they are both based on the PCBC (Fig. 1).While the model of Cuccia et al. is based on the PCBC, it does not overlap with the Hankel Transform of R PCBC ðρ; z 0 Þ.The derived model R PCBC ðkÞ does overlap with the Hankel Transform of R PCBC ðρ; z 0 Þ.An example of the reflectance versus μ 0 s f −1 for the different models and MC simulations is shown in Fig. 2 for μ a ¼ 0.01 mm −1 and μ 0 s ¼ 5 mm −1 .For higher spatial frequencies, measurements are in the subdiffuse regime, where the total reflectance should be equal to the sum of the semiballistic and diffuse reflectance. 18Thus, for any value of μ 0 s f −1 the diffuse reflectance should be equal to or smaller than the total reflectance.In Fig. 2, this is only true for R PCBC , which indicates that both R Cuccia and R EBC overestimate the diffuse reflectance.
Figure 3 depicts the median and interquartile ranges of the relative error in the reflectance for the three different models compared to all MC simulations.To better show the difference between the results for R Cuccia and R PCBC Fig. 3(b) shows a zoom-in view.Using R PCBC ðkÞ instead of R Cuccia ðkÞ reduces the median relative error by 45% from 0.011 to 0.006.R EBC ðkÞ results in a larger median relative error of 0.018, and also a much larger range of errors compared to the other two models.
While Fig. 3 compares the relative errors in the expected reflectance, we also determined the relative errors for each model in extracted values of μ a and μ s 0 from reflectance values obtained from MC simulations for spatial frequencies of 0 and of 0.5 mm −1 (Fig. 4).Using R PCBC ðkÞ instead of R Cuccia ðkÞ reduces the median relative error in the extracted value of μ a by 64% from 0.022 to 0.008 and the median relative error in the extracted value of μ 0 s by 10% from 0.029 to 0.026.The highest median relative errors were obtained with R EBC ðkÞ: 0.048 for μ a and 0.085 for μ s 0 .Also, the distribution of errors is the largest for R ECB ðkÞ.10), black lines] for two sets of optical properties as indicated in each subfigure.R Cuccia does not match the Hankel transform of Eq. ( 6) for the R PCBC ðρ; z 0 Þ, while R PCBC ðk Þ does match.
Fig. 2 (a) Models for the diffuse reflectance compared to the total simulated reflectance for μ a ¼ 0.01 mm −1 and μ 0 s ¼ 1 mm −1 for the HG phase function with g 1 ¼ 0.9 and (b) the two-term HG phase function with g 1 ¼ 0.83.In both, the diffuse (higher values of μ 0 s f −1 ) and the subdiffuse regime (lower values of μ 0 s f −1 ) the diffuse reflectance should be equal to or lower than the total reflectance.This is only true for R PCBC in this figure.Thus, R Cuccia and R EBC overestimate the diffuse reflectance.

Discussion
We developed a new model for the diffuse reflectance in SFDI that reduces the error in the estimated reflectance and the extracted optical properties compared to the currently used model of Cuccia et al. 2 Our model was derived by integrating the response function for a pencil illumination RðρÞ over ρ from 0 to ∞ under the PCBC for an extended source (yielding the response for a spatially unmodulated source) and replacing μ eff by μ 0 to take the influence of spatially modulated illumination into account.In other words, we hypothesized and demonstrated the equivalence between (a) computing the Hankel transform of a pencil beam response RðρÞ as in Eq. ( 5) to (b) direct integration of RðρÞ over ρ from 0 to ∞ followed by the substitution of μ eff by μ 0 eff .We compared the resulting equations for two models for RðρÞ based on two different boundary conditions, the PCBC and the EBC.We found that the errors in the predicted reflectance, as well as in extracted optical properties were much lower for the PCBC than the EBC.More importantly, for high spatial frequencies measurements are in the subdiffuse regime, where the total reflectance is the sum of a diffuse and a semiballistic component.Thus, for high spatial frequencies, the diffuse reflectance should always be lower than the total reflectance.For the EBC and the model of Cuccia et al., this was not the case.There is no theoretical reason to prefer our approach (starting from a pencil beam illumination) to the approach of Cuccia et al. (starting from a plane wave illumination), or to prefer the PCBC over the EBC as all are based on sound physical principles. 12The differences in the reflectance values obtained with the different approaches demonstrate the limitations in modeling diffuse light transport in general.Solving this apparent ambiguity is beyond the scope of this manuscript.
][21][22][23][24] In this paper, we used D ¼ 1∕ð3ðμ s 0 þ μ a ÞÞ to ensure a fair comparison to the model of Cuccia et al., which used absorption in the diffusion coefficient.In the diffuse regime, the reduced scattering coefficient is much larger than the absorption coefficient and the reflectance values obtained with each model would, thus, barely be different regardless of whether or not the absorption coefficient was included in the diffusion coefficient.Therefore, for the purpose of this paper, we could not determine which definition of the diffusion coefficient is more appropriate for SFDI.Even so, the definition of the diffusion coefficient does become important for the development of a subdiffuse model for SFDI.If the subdiffuse reflectance would be modeled as the sum of a diffuse and a semiballistic term, 18 the accuracy of the diffuse term for high values of μ a is important.Currently, a model does exist for subdiffuse SFDI, 7 but it is only valid for one type of tissue phase function and does not include absorption.Therefore, the currently available subdiffuse model cannot be used to interrogate tissue, since there will always be absorbers present and the tissue phase function is generally not known.
Apart from analytical models, approaches exist to extract optical properties from SFDI measurements based on MC simulations-either employing look-up-tables or machine-learning algorithms.MC simulations can include the details of the optical setup that is used (such as angle of incidence and detector numerical aperture) and look-up-tables and machine-learning algorithms can improve the speed of optical property determination.For medical applications where speed is essential, such as endoscopy, 25 fast algorithms or look-up-tables to extract optical properties are favorable.However, when speed is less important, analytical models provide a few benefits.First, while machine-learning algorithms could overfit the problem and might not provide accurate answers for optical properties that were not simulated, this is much less likely for analytical models.Second, while anybody can use an analytical model since it is a formula that is written out, machine-learning algorithms are often not freely available and can only be reproduced by redoing the MC simulations and retraining the algorithm.Regardless, analytical models are valuable for our fundamental understanding of light transport in general and for SFDI specifically.For example, previously we developed an analytical model for another subdiffuse spectroscopy technique that uses fiber-optic probes: single fiber reflectance (SFR) spectroscopy.We modeled the subdiffuse reflectance as the sum of a diffuse and a semi-ballistic component and identified the new parameter p sb to incorporate the influence of the phase function on the semi-ballistic component. 26Without the analytical model that we developed for the diffuse reflectance in SFR spectroscopy. 27we would not have been able to identify the parameter p sb .More importantly, now that we determined that the subdiffuse reflectance in SFR spectroscopy depends on μ 0 s , μ a , and p sb , it is possible to properly perform MC simulations to create look-uptables or machine-learning algorithms for SFR spectroscopy by including simulations with a range of μ 0 s , μ a and p sb values.In conclusion, we investigated the diffuse reflectance in SFDI by comparing the currently used model of Cuccia et al. to two new models based on integrating the theoretical response function RðρÞ under the EBC and PCBC for a pencil beam illumination over ρ from 0 to ∞ and replacing μ eff by μ eff 0 .The model based on the PCBC provides the best results, and reduces the median relative error by 10% for the extracted μ 0 s , 64% for μ a and 45% for the reflectance.Errors in the expected reflectance can further influence the accuracy of extracted optical properties, since SFDI measurements involve a calibration procedure with a phantom with known optical properties for which expected reflectance is used to calibrate the setup.

Fig. 1
Fig. 1 (a) and (c) Reflectance versus μ 0 s f −1 calculated with the Hankel transform of Eq. (6) for the R PCBC ðρ; z 0 Þ (red lines) compared to R Cuccia [Eq.(4), blue lines] and (b) and (d) the model for R PCBC ðk Þ [Eq.(10), black lines] for two sets of optical properties as indicated in each subfigure.R Cuccia does not match the Hankel transform of Eq. (6) for the R PCBC ðρ; z 0 Þ, while R PCBC ðk Þ does match.

Fig. 4
Fig. 4 Median and interquartile ranges of the relative error in the extracted values for μ a and μ 0 s with R Cuccia ðk Þ [Eq.(4)], R PCBC ðk Þ [Eq.(10)], and R EBC ðk Þ [Eq.(11)].Using R PCBC ðk Þ instead of R Cuccia ðk Þ reduces the median relative error in the extracted value of μ a by 64% from 0.022 to 0.008 and the median relative error in the extracted value of μ 0 s by 10% from 0.029 to 0.026.