1.

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,³ such as tumor margin assessment in breast cancer⁴ and assessment of diabetic feet.⁵ In 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_{\mathrm{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⁶ to match the clinical application.

To extract optical properties from SFDI measurements, the reflected intensity is first demodulated. The amplitude modulation $M_{\mathrm{AC}}$ for any tissue location $x$, and projected spatial frequency $f_x$ can be expressed as

To obtain $R_{\text{tissue}}$ from $M_{\mathrm{AC}}$, the modulation transfer function of the optical system has to be determined. ${\mathrm{MTF}}_{\mathrm{SYS}}$ is obtained through a calibration step where the $M_{\mathrm{AC}}$ of a reference sample is measured for which the optical expected reflectance $R_{\mathrm{ref}}$ is known,

To extract tissue optical properties from SFDI measurements of $R_{\text{tissue}}$, two main approaches exist. The first approach is to use an analytical model based on physics that relates $R_{\text{tissue}}$ to the tissue optical properties.^2,7 The second approach is to use Monte Carlo (MC) simulations to relate $R_{\text{tissue}}$ to tissue optical properties, either through a look-up-table approach or based on machine learning.^8–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² and the model of Kanick et al. for the subdiffuse regime.⁷ The 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\ll {\mu }_{\mathrm{tr}}$ and ${\mu }_s^{\prime }\gg {\mu }_a$ (where ${\mu }_{\mathrm{tr}}={\mu }_a+{\mu }_s^{\prime }$). 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.

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(\rho )$, 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:

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(\rho )$ 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(\rho )$ 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.¹² and the extended boundary condition (EBC) as proposed by Farrell et al.¹³ For the PCBC, the irradiance at the boundary is set equal to the integral of the reflected radiance¹² 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.¹⁴ The two boundary conditions result in the following expressions for the reflectance:

The reflectance as a function of radial distance is then found as $R(\rho )={\int }_{z_0=0}^{\infty }R(\rho ,z_0)S(z_0){\mathrm{dz}}_0$, where $S$ is the source term, for which we use a distributed source term $S(z)=a^{\prime }{\mu }_{\mathrm{tr}}\exp (-{\mu }_{\mathrm{tr}}z)$. To obtain our new models for the SFDI reflectance we integrate $R(\rho )$ for each boundary condition over the radial coordinate $\rho$ from 0 to $\infty$ and we replace ${\mu }_{\mathrm{eff}}$ with ${\mu }_{\mathrm{eff}}^{\prime }$.

Performing the integrations for $R_{\mathrm{PCBC}}(\rho ,z_0)$ and $R_{\mathrm{EBC}}(\rho ,z_0)$ over $z_0$ and $\rho$ we arrive at

3.

Methods

Since the model of Cuccia et al. was based on the PCBC, we first compare $R_{\text{Cuccia}}(k)$ and $R_{\mathrm{PCBC}}(k)$ to the Hankel Transform of $R_{\mathrm{PCBC}}(\rho ,z_0)$. To calculate this Hankel Transform with Matlab we first integrated $R_{\mathrm{PCBC}}(\rho ,z_0)$ over ${\mathrm{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_{\text{Cuccia}}(k)$, $R_{\mathrm{PCBC}}(k)$ and $R_{\mathrm{EBC}}(k)$) we performed MC simulations to obtain the reflectance versus radial distance, $R_{\mathrm{MC}}(\rho )$. To calculate the reflectance measured by SFDI, $R_{\mathrm{MC}}(k)$, we performed a Hankel Transform. For each model, we calculated the relative error (RE) in the reflectance with respect to $R_{\mathrm{MC}}(k)$ for all the simulations

We also determined the relative errors in the values of ${\mu }_{\mathrm{a}}$ and ${{\mu }_{\mathrm{s}}}^{\prime }$ 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{\mathrm{mm}}^{-1}$.

4.

Results

We first compared $R_{\text{Cuccia}}(k)$ [Eq. (4)] and $R_{\mathrm{PCBC}}(k)$ [Eq. (10)] to the Hankel Transform [Eq. (5)] of $R_{\mathrm{PCBC}}(\rho ,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_{\mathrm{PCBC}}(\rho ,z_0)$. The derived model $R_{\mathrm{PCBC}}(k)$ does overlap with the Hankel Transform of $R_{\mathrm{PCBC}}(\rho ,{\mathrm{z}}_0)$.

Fig. 1

(a) and (c) Reflectance versus ${\mu }_s^{\prime }{f}^{-1}$ calculated with the Hankel transform of Eq. (6) for the $R_{\mathrm{PCBC}}(\rho ,z_0)$ (red lines) compared to $R_{\text{Cuccia}}$ [Eq. (4), blue lines] and (b) and (d) the model for $R_{\mathrm{PCBC}}(k)$ [Eq. (10), black lines] for two sets of optical properties as indicated in each subfigure. $R_{\text{Cuccia}}$ does not match the Hankel transform of Eq. (6) for the $R_{\mathrm{PCBC}}(\rho ,z_0)$, while $R_{\mathrm{PCBC}}(k)$ does match.

An example of the reflectance versus ${\mu }_s^{\prime }{f}^{-1}$ for the different models and MC simulations is shown in Fig. 2 for ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=5{\mathrm{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.¹⁸ Thus, for any value of ${\mu }_s^{\prime }{f}^{-1}$ the diffuse reflectance should be equal to or smaller than the total reflectance. In Fig. 2, this is only true for $R_{\mathrm{PCBC}}$, which indicates that both $R_{\text{Cuccia}}$ and $R_{\mathrm{EBC}}$ overestimate the diffuse reflectance.

Fig. 2

(a) Models for the diffuse reflectance compared to the total simulated reflectance for ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ for the HG phase function with ${\mathrm{g}}_1=0.9$ and (b) the two-term HG phase function with ${\mathrm{g}}_1=0.83$. In both, the diffuse (higher values of ${\mu }_{\mathrm{s}}^{\prime }{f}^{-1}$) and the subdiffuse regime (lower values of ${\mu }_{\mathrm{s}}^{\prime }{f}^{-1}$) the diffuse reflectance should be equal to or lower than the total reflectance. This is only true for $R_{\mathrm{PCBC}}$ in this figure. Thus, $R_{\text{Cuccia}}$ and $R_{\mathrm{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_{\text{Cuccia}}$ and $R_{\mathrm{PCBC}}$ Fig. 3(b) shows a zoom-in view. Using $R_{\mathrm{PCBC}}(k)$ instead of $R_{\text{Cuccia}}(k)$ reduces the median relative error by 45% from 0.011 to 0.006. $R_{\mathrm{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.

Fig. 3

(a) The median and interquartile ranges of the relative error in $R_{\text{Cuccia}}(k)$ [Eq. (4)], $R_{\mathrm{PCBC}}(\mathrm{k})$ [Eq. (10)] and $R_{\mathrm{EBC}}(k)$ [Eq. (11)] versus the simulated reflectance. (b) Zoom-in view of Fig. 3(a). Using $R_{\mathrm{PCBC}}(k)$ instead of $R_{\text{Cuccia}}(k)$ reduces the median relative error in the expected reflectance by 45% from 0.011 to 0.006.

While Fig. 3 compares the relative errors in the expected reflectance, we also determined the relative errors for each model in extracted values of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ from reflectance values obtained from MC simulations for spatial frequencies of 0 and of $0.5{\mathrm{mm}}^{-1}$ (Fig. 4). Using $R_{\mathrm{PCBC}}(k)$ instead of $R_{\text{Cuccia}}(k)$ reduces the median relative error in the extracted value of ${\mu }_a$ by 64% from 0.022 to 0.008 and the median relative error in the extracted value of ${\mu }_s^{\prime }$ by 10% from 0.029 to 0.026. The highest median relative errors were obtained with $R_{\mathrm{EBC}}(\mathrm{k})$: 0.048 for ${\mu }_a$ and 0.085 for ${{\mu }_s}^{\prime }$. Also, the distribution of errors is the largest for $R_{\mathrm{ECB}}(k)$.

Fig. 4

Median and interquartile ranges of the relative error in the extracted values for ${\mu }_a$ and ${\mu }_s^{\prime }$ with $R_{\text{Cuccia}}(k)$ [Eq. (4)], $R_{\mathrm{PCBC}}(k)$ [Eq. (10)], and $R_{\mathrm{EBC}}(k)$ [Eq. (11)]. Using $R_{\mathrm{PCBC}}(k)$ instead of $R_{\text{Cuccia}}(k)$ reduces the median relative error in the extracted value of ${\mu }_a$ by 64% from 0.022 to 0.008 and the median relative error in the extracted value of ${\mu }_{\mathrm{s}}^{\prime }$ by 10% from 0.029 to 0.026.

5.

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.² Our model was derived by integrating the response function for a pencil illumination $R(\rho )$ over $\rho$ from 0 to $\infty$ under the PCBC for an extended source (yielding the response for a spatially unmodulated source) and replacing ${\mu }_{\mathrm{eff}}$ by ${\mu }_{\mathrm{eff}}^{\prime }=\sqrt{{\mu }_{\mathrm{eff}}^2+k^2}$ 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(\rho )$ as in Eq. (5) to (b) direct integration of $R(\rho )$ over $\rho$ from 0 to $\infty$ followed by the substitution of ${\mu }_{\mathrm{eff}}$ by ${\mu }_{\mathrm{eff}}^{\prime }$.

We compared the resulting equations for two models for $R(\rho )$ 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.¹² The 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.

There has been a debate in literature on whether or not the absorption coefficient should be incorporated in the diffusion coefficient for the analytic solutions to the PCBC and EBC.^19–24 In this paper, we used $D=1/(3({{\mu }_s}^{\prime }+{\mu }_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,¹⁸ the accuracy of the diffuse term for high values of ${\mu }_{\mathrm{a}}$ is important. Currently, a model does exist for subdiffuse SFDI,⁷ 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,²⁵ 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_{\mathrm{sb}}$ to incorporate the influence of the phase function on the semi-ballistic component.²⁶ Without the analytical model that we developed for the diffuse reflectance in SFR spectroscopy.²⁷ we would not have been able to identify the parameter $p_{\mathrm{sb}}$. More importantly, now that we determined that the subdiffuse reflectance in SFR spectroscopy depends on ${\mu }_s^{\prime }$, ${\mu }_a$, and $p_{\mathrm{sb}}$, it is possible to properly perform MC simulations to create look-up-tables or machine-learning algorithms for SFR spectroscopy by including simulations with a range of ${\mu }_s^{\prime }$, ${\mu }_a$ and $p_{\mathrm{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(\rho )$ under the EBC and PCBC for a pencil beam illumination over $\rho$ from 0 to $\infty$ and replacing ${\mu }_{\mathrm{eff}}$ by ${{\mu }_{\mathrm{eff}}}^{\prime }$. The model based on the PCBC provides the best results, and reduces the median relative error by 10% for the extracted ${\mu }_s^{\prime }$, 64% for ${\mu }_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 the expected reflectance is used to calibrate the setup.