1.

Introduction

The skin is one of the largest organs of the human body, and it acts as an interface to the outside world and as a barrier to protect from external pathogens and other hazards. A wound is formed when skin is damaged or otherwise unable to perform its functions.^1,2 The most common method used for evaluating the progression of the healing process is visual inspection by a trained physician. Parameters, such as wound size, color, granulation, and presence of secretions or necrotic tissue are annotated at different points in time and used to make an estimate of how much a wound has progressed.^3–5 This method is fast, repeatable, and non-invasive, but the dependence on the expertise of the clinical personnel makes it a non-reliable way to obtain a consistent diagnosis.⁶ A more specific method to obtain information on the wound is biopsy, followed by staining and cell histology. This approach gives more structural information at the cellular level, but it is invasive, extremely localized and non-repeatable (especially on patients in-vivo), so it is only performed when there is already suspicion of malignant tissue. An ideal technique for wound assessment would combine the advantages of both visual inspection and histology, being non-invasive, repeatable in time, and specific, while presenting an objective measurement of parameters that can be employed as a metric for the evaluation of the efficiency of different regeneration treatments.

Optical techniques are non-invasive, fast methods for performing measurements on tissue using light. They satisfy most of the prerequisites that we want from an ideal wound healing assessment method: from the way that light interacts with tissue, it is possible to quantitatively measure concentrations of biological chromophores and indirectly obtain structural information on a microscopic scale, over a large field of view. Optical techniques are already widely utilized in clinical practice, where they are used for determining blood perfusion and blood oxygen saturation,^7–9 detecting cancer,^10,11 and optimizing photodynamic therapy,^12,13 among other applications.^14,15

Spatial frequency domain imaging (SFDI) is an optical imaging technique that makes use of structured light illumination to quantitatively measure the absorption (${\mu }_a$) and reduced scattering (${\mu }_s^{\prime }$) coefficients in tissue and separate their effects by making use of the frequency-specific response of the tissue.¹⁶ The absorption of light contains information about the concentration of molecules that are capable of absorbing photons (e.g., melanin, haemoglobin, etc.), whereas light scattering is dependent on the density and size distribution of the microscopic structures (e.g., collagen fibers, cell nuclei, etc.). These properties make SFDI a potentially valuable tool for the diagnosis of wound healing as it is able to measure objective parameters related to tissue function (e.g., blood concentration and oxygenation) and morphology (e.g., cellular proliferation and tissue remodeling) in a fast, non-invasive, and repeatable way. The technique, in its current state, has a few limitations, most notable of which is that the models of light transport, used to quantify the absorption and scattering parameters, presume a semi-infinite, optically homogeneous geometry. When looking at the practical application of wound assessment, it is evident that a homogeneous model is not an accurate representation of the underlaying physiology: biological tissue is quite heterogeneous and contains multiple thin layers with different properties. For the purpose of overcoming the limitations of homogeneous models we introduce a two-layer model composed of a thin layer, with a thickness in the range of hundreds of microns (which simulates the new tissue growth) and a thick layer, which can be approximated by a semi-infinite geometry and simulates the underlaying wound site. The introduction of the two-layer model alone, however, is not sufficient. To discriminate between the optical properties between the top and the bottom layer, we also need a dataset that contains information spanning multiple penetration depths of light into the tissue volume. It was shown in previous investigations that, when acquiring data over wavelengths of light spanning the visible and near infrared, each will have a different penetration depth depending on the optical properties of the tissue, with longer wavelengths (e.g., near-infrared) normally having longer penetration depths compared with shorter wavelengths (e.g., visible).¹⁷

However, in the context of wound healing, this approach is not as applicable because the absorption contrast between layers is less pronounced. Scattering contrast, on the other hand, would contain distinct information about depth-specific structural changes in the wound site. SFDI also has a unique property of allowing for varying the penetration depth of the light patterns independently from the wavelength of light. This is made possible by means of the spatial frequency ($f_x$) of the illumination patterns.¹⁸ Sinusoidal patterns with higher frequencies have a lower penetration depth compared with patterns with lower frequencies. By acquiring multiple datasets with an increasing average spatial frequency, we are effectively sampling smaller and smaller volumes, which contain information in increasing proportion from the superficial layer, compared with the bottom layer. This multi-frequency approach is particularly suited to detecting variations in scattering contrast over thin superficial layer, as absorption contrast is discernible mostly at low $f_x$ (and higher penetration depth), whereas scattering has a larger influence at high $f_x$ (smaller penetration depth).

In this study, we acquire an SFDI dataset at several spatial frequency and then subdivide it in multiple overlapping sub-sets containing increasing spatial frequencies. These sub-sets, containing data with different penetration depth by means of different $f_x$, constitute the basis of the two-layered model proposed in this study and are used to determine the layer specific scattering coefficients, which for the practical application of wound assessment can be correlated to new cell growth (re-epithelialization) and morphological changes within the underlying wound bed (remodeling). To test the validity of the model, SFDI data were acquired on silicone multi-layered optical phantoms with known properties and then processed according to the previously described multi-frequency approach. These tissue simulating multi-layer phantoms were designed to mimic the relative, layer-specific optical properties observed during re-epithelialization of a wound. Using the individual properties of the two layers of the phantoms and their geometrical parameters, we compared the experimental data with the proposed two-layer model using three different analytical models of light transport. The root mean square percentage error (RMSPE) of the modeled data was then calculated as an evaluation metric to determine the performance of the models in different ranges of parameters.

2.

Two-layer Scattering Model

Thin layered models of skin, based on the difference between the absorption in the epidermis and dermis layers, are already used in clinical applications to account for the different absorption levels of epidermis and dermis.^17,19 Scattering-based multi-layer models, however, are less common. SFDI can be used to measure both absorption and scattering in a tissue, but for the purpose of this study, we consider differences in the scattering properties between the layers. As mentioned, in the context of wound healing, the difference in absorption between the layers is less pronounced, whereas the scattering coefficient can be subjected to large variations that give insight on the changes in tissue morphology due to the healing process. A two-layer scattering model based on the spectral difference in depth penetrance between the layers was already developed by our group to investigate morphologic changes within skin pigmentation and melanoma; however, that approach was optimized to isolate tissue layer thicknesses of $\sim 0.2$ to 4 mm.²⁰ In the present work, we aim to mimic the physiology of a healing wound, in which sub-cellular organelles of different sizes and new epithelial cells give different sources of scattering-based contrast. The end goal is to be able to measure this scattering contrast between thin tissue layers (0.05 to 0.5 mm).

The geometry of our two-layer model is represented in Fig. 1. The thin top layer has a finite thickness d and a scattering coefficient ${\mu }_{S_{\text{top}}}^{\prime }$. The bottom layer is a semi-infinite slab, with scattering coefficient ${\mu }_{S_{\mathrm{bot}}}^{\prime }$. The index of refraction ($n$) is also assumed to be the same in both layers to avoid unwanted reflections and refractions due to index mismatch.

Fig. 1

2D drawing of the geometry of the two-layer scattering model (left). A thin layer of thickness $d$ with scattering coefficient ${\mu }_{S_{\mathrm{top}}}^{\prime }$ is laid on a semi-infinite layer with scattering coefficient ${\mu }_{S_{\mathrm{bot}}}^{\prime }$. Both layers are assumed to be homogeneous and extend infinitely in the $xy$ plane. The light fluence $\varphi$ (right) decreases exponentially from the surface in the $z$ direction (depth). The weighting coefficient $\alpha (d)$ is defined as the partial contribution of fluence squared over the thin layer of thickness $d$. The contribution is obtained by integrating ${\varphi }^2$ up to $d$ and normalizing to the total integral of ${\varphi }^2$.

There are three parameters in this scattering model: ${\mu }_{S_{\text{top}}}^{\prime }$, ${\mu }_{S_{\mathrm{bot}}}^{\prime }$, and $d$. When performing normal SFDI measurements, we obtain a single measured scattering coefficient (${\mu }_{s_{\mathrm{meas}}}^{\prime }$), which is assumed to derive from a single, homogeneous, and semi-infinite geometry. We want to define a modeled scattering coefficient (${\mu }_{s_{\mathrm{mod}}}^{\prime }$) that models the behavior of ${\mu }_{s_{\mathrm{meas}}}^{\prime }$ and takes into consideration the relative contributions of the three parameters of the two-layer model. A first assumption that we make is that ${\mu }_{s_{\mathrm{mod}}}^{\prime }$ is a linear combination of the two free scattering parameters, ${\mu }_{S_{\text{top}}}^{\prime }$ and ${\mu }_{S_{\mathrm{bot}}}^{\prime }$, given as

2.3.

Analytical Models of Fluence ($\varphi$)

The last part needed to define our two-layer model is to choose a suitable model of light fluence ($\varphi$) that solves the radiative transport equation, specific to the spatial frequency domain. There are several options available to calculate a statistically reliable representation of $\varphi$: in general, from analytical models that make use of approximations to obtain simplified mathematical formulations of $\varphi$^23–25 to computer simulations using the Monte Carlo method that use realistic geometries and random propagation for millions of photons.^26,27 In this study, we selected and evaluated the performance of two analytical models of light fluence. This decision was made because of the ease of implementation and much faster computation times compared with simulation-based models. In addition, normal Monte Carlo methods cannot currently perform direct simulations in the spatial frequency domain (SFD), so the spatial frequency dependency is obtained indirectly by means of the Hankel transform,²⁸ which is only applicable in certain conditions, at the cost of losing all spatial information.

2.3.1.

Standard diffusion approximation

The first analytical model that we investigated is the standard diffusion approximation (SDA).²³ The SDA is a model already used in many medical applications that works best in regimes in which scattering is predominant compared with absorption (${\mu }_s^{\prime }\gg {\mu }_a^{\prime }$).¹⁶ For this reason, the SDA is more effective for light in the infrared spectrum because of the low absorption from biological tissue. The formula of light fluence (normalized to the incident power $P_0$) in the SFD for the diffuse approximation is given in the following equation:

Eq. (4)

$$\frac{\varphi (z,f_x)}{P_0}=A·\exp (-{\mu }_{tr}·z)+C·\exp (-{\mu }_{\mathrm{eff}}^{\prime }(f_x)·z),$$

where $A$ and $C$ are constants derived by the choice of an appropriate boundary condition and $R$ is a constant dependent on the effective reflection coefficient, which are given as

Eq. (5)

$$A=\frac{3\frac{{\mu }_s^{\prime }}{{\mu }_{tr}}}{\frac{{\mu }_{eff}^{\prime 2}}{{\mu }_{tr}^2}-1};C=\frac{-3\frac{{\mu }_s^{\prime }}{{\mu }_{tr}}(1+3R)}{(\frac{{\mu }_{\mathrm{eff}}^{\prime 2}(f_x)}{{\mu }_{tr}^2}-1)(\frac{{\mu }_{\mathrm{eff}}^{\prime }(f_x)}{{\mu }_{tr}}+3\mathrm{R})};R=\frac{1-R_{\mathrm{eff}}}{2(1+R_{\mathrm{eff}})}.$$

For the full derivation and explanation of the coefficients, we refer to Cuccia et al.¹⁶

2.3.2.

$\delta$-P1 approximation

Because our SFDI instrumentation uses a light source that includes the visible spectrum and we are trying to detect photons that have short pathlengths, the SDA becomes less accurate as we are approaching the limits at which the model is still reliable. For this reason, we also considered a second analytical model. The $\delta$-P1 approximation is a diffusion-based model that introduces a correction factor to account for photons that are non-diffuse. This correction factor allows for modeling fluence more accurately in those conditions in which the SDA becomes less reliable (e.g., for short distance photon propagation).^24,29 The original $\delta$-P1 model did not include dependency on the spatial frequency $f_x$,²⁹ which was only introduced later, in the context of a doctoral thesis.³⁰ The formula of light fluence (normalized to the incident power $P_0$) in the spatial frequency domain for the $\delta$-P1 model is given in Eq. (6). For a complete derivation and explanation of the coefficients, we refer to the doctoral thesis of Seo.³⁰

Eq. (6)

$$\frac{{\varphi }_{AC}(z,f_x)}{P_0}=\frac{C^{*}}{{\mu }_{\mathrm{eff}}^{\prime }-{\mu }_{tr}^{*}}[\exp (-{\mu }_{tr}^{*}z)-\exp (-{\mu }_{\mathrm{eff}}^{\prime }z)]+\frac{C^{*}}{{\mu }_{\mathrm{eff}}^{\prime }+{\mu }_{tr}^{*}}[\exp (-{\mu }_{tr}^{*}z)-\exp (-{\mu }_{\mathrm{eff}}^{\prime }(z+2z_b))],$$

where ${\mu }_{tr}^{*}={\mu }_a+{\mu }_s^{*}={\mu }_a+(1-g^2){\mu }_s$ as described in Carp et al.²⁹ and the constants $C^{*}$ and $z_b$ are obtained by the application of the boundary condition, as described in the appendix of Seo’s doctoral thesis and are given as³⁰

Eq. (7)

$$C^{*}=\frac{3{\mu }_{tr}{\mu }_s^{*}}{2{\mu }_{\mathrm{eff}}^{\prime }};z_b=\frac{2}{3{\mu }_{tr}}\frac{(1+R_1)}{(1-R_1)}.$$

2.3.3.

Modified $\delta$-P1 approximation

In this study, we also derived our own spatial frequency-dependent $\delta \text{-}\mathrm{P}1$ equation based on the original $\delta \text{-}\mathrm{P}1$ equation²⁹ and applied the same modification and boundary conditions made by Cuccia et al. in their derivation of the SDA in the spatial frequency domain.¹⁶ We call this derivation modified $\delta \text{-}\mathrm{P}1$, or mod-$\delta \text{-}\mathrm{P}1$. The equation of light fluence for the mod-$\delta \text{-}\mathrm{P}1$ approximation is shown in the following equation:

Eq. (8)

$$\frac{\varphi (z,f_x)}{P_0}=(1+A^{\prime })\exp (-{\mu }_{tr}^{*}·z)+C^{\prime }\exp (-{\mu }_{\mathrm{eff}}^{\prime }z),$$

where the coefficients $A^{\prime }$ and $C^{\prime }$ are obtained by solving Boltzmann differential equation with the appropriate boundary conditions as

Eq. (9)

$$A^{\prime }=\frac{3{\mu }_s^{*}({\mu }_{tr}^{*}+g^{*}{\mu }_a)}{({\mu }_{\mathrm{eff}}^{\prime 2}-{\mu }_{tr}^{*2})};C^{\prime }=\frac{-A^{\prime }(1+\frac{2}{3}{R}^{\prime }{\mu }_{tr}{\mu }_{tr}^{*})+2R^{\prime }{\mu }_{tr}{g}^{*}{\mu }_s^{*}}{(1+\frac{2}{3}{R}^{\prime }{\mu }_{tr}{\mu }_{\mathrm{eff}}^{\prime })};R^{\prime }=\frac{(1+R_1)}{(1-R_1)}.$$

For a full derivation and explanation of the coefficients, we refer to Appendix A.

A comparison between the SDA and the two $\delta \text{-}\mathrm{P}1$ models of fluence are given in Fig. 2, which shows the values of $\varphi (z)$, normalized to the total integral of $\varphi$ over $z$. The models were calculated for optical properties ${\mu }_a=0.05{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=5{\mathrm{mm}}^{-1}$. Figure 2 also shows how the spatial frequency changes the distribution of the photons, with higher values of $f_x$ having photons more concentrated near the surface.

Fig. 2

Comparison between fluence ($\varphi$) models using SDA (left), the original $\delta \text{-}\mathrm{P}1$ approximation (center), and our mod-$\delta$-P1 derivation (right). The values of $\varphi$ were normalized to their total integral over $z$. In each plot, the effect of the spatial frequency ($f_x$) on the fluence distribution is also shown. All models of $\varphi$ in the figure were simulated on an homogeneous geometry with the same optical properties: ${\mu }_a=0.05{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=5{\mathrm{mm}}^{-1}$.

3.

Materials and Methods

3.3.

Data Processing

Our new multi-frequencies processing approach consists of subdividing the dataset into several smaller, partially overlapping sub-sets containing four frequencies each, with increasing values of $f_x$, as seen in Fig. 3. Each sub-set is individually processed to obtain ${\mu }_a$ and ${\mu }_s^{\prime }$ values, which have different penetration depths, as defined in Sec. 2. The average $f_x$ of each subset is used in the fluence calculation, denoted as $⟨f_x⟩$ in the text. The SFDI data were processed assuming a homogeneous volume using the procedure described in Cuccia et al.:¹⁶ the amplitude of the modulated patterns was extracted and a homogeneous reference phantom with known optical properties was used for calibration to obtain the measurements of diffuse reflectance ($R_d$). A white Monte Carlo model of $R_d$ is then used in an iterative algorithm that changes the input parameters $({\mu }_a,{\mu }_s^{\prime })$ according to an optimization strategy, until the $R_d$ of the model matches the $R_d$ measured on the target, giving as a result the $({\mu }_a,{\mu }_s^{\prime })$ pair of the measured target.

Fig. 3

Subdivision of the dataset containing 11 spatial frequencies ($f_x$) in 8 smaller, partially overlapped sub-sets, containing $4f_x$ each. The sub-sets have different penetration depths ($\delta$), which are estimated considering the average $f_x$ of the sub-sets ($⟨f_x⟩$), with lower $⟨f_x⟩$ resulting in higher $\delta$.

4.

Results

4.1.

Measured ${\mu }_s^{\prime }$

A sample of the data measured in Sec. 3.2.2 and processed assuming homogeneous tissue, as described in Sec. 3.3, is shown in Fig. 4. The figure contains data from a single spectral band (536 nm) because the measurements at each wavelength are independent from one another. The graph shows the scattering coefficients of the five two-layered phantoms and how they change with $⟨f_x⟩$. The horizontal axis is the average of the spatial frequencies contained in each of the eight sub-sets. The ${\mu }_s^{\prime }$ of the top layer (${\mathrm{TiO}}_2$ – light orange line) and bottom layer (${\mathrm{Al}}_2{\mathrm{O}}_3$ – dark blue line) are included for comparison. The ${\mu }_s^{\prime }$ of these two layers are assumed to be constant across the spatial frequencies and were obtained by measuring it on the thick homogeneous phantoms using the same SFDI system.

Fig. 4

Scattering coefficients (${\mu }_s^{\prime }$) of the five two-layered phantoms (colored markers), measured at 536 nm. The solid lines are the ${\mu }_s^{\prime }$ of the top layer (light orange) and bottom layer (dark blue). The data on the horizontal axis represent the average spatial frequency $⟨f_x⟩$ of the eight sub-sets into which the data were divided.

4.2.

Fluence-Based Models

Using the values of ${\mu }_{s_{\text{top}}}^{\prime }$ and ${\mu }_{s_{\mathrm{bot}}}^{\prime }$ measured in Sec. 3.2.1 and the phantom thickness ($d$) reported in Table 1, the expected values of ${\mu }_s^{\prime }(f_x)$ of the two-layered phantoms were modelled using Eqs. (1) and (3). The ${\mu }_a$ and ${\mu }_s^{\prime }$ measured on the layered phantoms, as described in Sec. 3.2.2, were used to calculate the fluence $\varphi$ used in Eq. (3) with the three models presented in Sec. 2.3. The procedure is summarized in the flowchart in Fig. 5. In Fig. 6, a comparison between the modeled ${\mu }_s^{\prime }$ values and the experimental measurement on phantoms at 536 nm is shown for both SDA and the two $\delta$-P1 models.

Fig. 5

Summary of the procedure to acquire depth-sensitive data using multi-frequency SFDI and use it as input to the two-layer scattering model. SFDI images are acquired at several spatial frequencies and divided in sub-sets containing four frequencies each. From each sub-set, optical properties are estimated and used to calculate the homogeneous fluence rate $\varphi$. Values of $\varphi$ are then used in combination with the model parameters to estimate the partial volume contributions $\alpha$ and obtain a model of ${\mu }_s^{\prime }(⟨f_x⟩)$.

Fig. 6

Comparison of the three fluence-based two-layer models (dashed lines) with the data (markers) measured at 536 nm. The top (light orange) and bottom (dark blue) solid lines are the scattering coefficients of the ${\mathrm{TiO}}_2$ and ${\mathrm{Al}}_2{\mathrm{O}}_3$ layers, respectively. The three models are (a) SDA, (b) the original $\delta \text{-}\mathrm{P}1$ approximation, and (c) our mod-$\delta \text{-}\mathrm{P}1$ approximation.

4.3.

Models’ Performance

To obtain an objective evaluation of the accuracy of the two-layered models of ${\mu }_s^{\prime }$ and compare their performance, the RMSPE was calculated with respect to the measurements described in Sec. 4.1 for each of the models and each of the phantom thicknesses using the following equation:

Eq. (10)

$$\mathrm{RMSPE}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(\frac{{\mu }_{s_{\mathrm{meas}}}^{\prime }(⟨f_i⟩)-{\mu }_{s_{\mathrm{mod}}}^{\prime }(⟨f_i⟩)}{{\mu }_{s_{\mathrm{meas}}}^{\prime }(⟨f_i⟩)}\right)}^2}*100,$$

where $N$ is the number of average spatial frequencies contained in the dataset ($N=8$), ${\mu }_{s_{\mathrm{meas}}}^{\prime }$ is the measured scattering coefficient, and ${\mu }_{s_{\mathrm{mod}}}^{\prime }$ is the scattering coefficient modeled using Eq. (1). The chart in Fig. 7 summarizes the RMSPE of three datasets having different scattering contrasts, defined as the ratio ${\mu }_{s_{\text{top}}}^{\prime }/{\mu }_{s_{\mathrm{bot}}}^{\prime }$.

Fig. 7

RMSPE calculated on the difference between the scattering coefficient of the three fluence-based two-layer models and the experimental measurements. Each figure has a different scattering contrast, defined as the ratio ${\mu }_{s_{\text{top}}}^{\prime }/{\mu }_{s_{\mathrm{bot}}}^{\prime }$ and shows the RMSPE of the three models for the five top layer thicknesses (in mm). The scattering ratios are approximately $a=4$, $b=3.5$, and $c=3$.

5.

Discussion

By observing the ${\mu }_{s_{\mathrm{meas}}}^{\prime }$ data in Fig. 4, we can make a few observations. First, the value of the ${\mu }_s^{\prime }$ coefficients measured on the two-layered phantoms are contained between ${\mu }_{s_{\text{top}}}^{\prime }$ and ${\mu }_{s_{\mathrm{bot}}}^{\prime }$, which is consistent with the hypothesis made in Sec. 2 of ${\mu }_{s_{\mathrm{mod}}}^{\prime }$ being a combination of the other two scattering coefficients. Second, with an increasing thickness of the thin layer, the ${\mu }_{s_{\mathrm{meas}}}^{\prime }$ of the two-layered phantoms becomes closer to ${\mu }_{s_{\text{top}}}^{\prime }$. This is to be expected because the relative contribution of the top layer to the overall scattering becomes greater the thicker the layer is. Similarly, with the increase of the average $f_x$, the ${\mu }_{s_{\mathrm{meas}}}^{\prime }$ of the two-layered phantoms becomes closer to ${\mu }_{s_{\text{top}}}^{\prime }$. This behavior is also in line with the SFDI property shown in Sec. 2.2, where higher spatial frequencies have less penetration depth, so the relative contribution from the top layer becomes higher. A last consideration that can be done is that the dependency of the two-layered ${\mu }_s^{\prime }$ over $⟨f_x⟩$ seems to be non-linear in nature and tends to “saturate” for high enough values of $⟨f_x⟩$ or layer thickness. This saturation of ${\mu }_s^{\prime }$ can be explained by assuming that the penetration depth of light ($\delta$) becomes smaller than the thickness of the top layer, so we are effectively measuring a single layer instead of two. Looking at the two-layer scattering models in Fig. 6, it is possible to determine where the underlying fluence models are the most accurate and what their limitations are. The SDA model [Fig. 6(a)] seems to be more accurate for larger thicknesses, but it tends to underestimate the contribution of the top layer, especially at low $⟨f_x⟩$. This is expected, as SDA assumes that the photons have long enough pathlengths to be considered diffuse (distance from the light source $\gg 1/{\mu }_t$). This behaviour is also evident from the fluence plots in Fig. 2 that show a larger concentration of photons deeper in the tissue compared to the other models. The original $\delta \text{-}\mathrm{P}1$ model [Fig. 6(b)] seems to have an accuracy similar to the SDA in the extreme cases ($d<0.1\mathrm{mm}$, $d>1\mathrm{mm}$), but it is more accurate in the middle range ($0.3\mathrm{mm}<d<0.7\mathrm{mm}$) and does not overestimate the contribution of the bottom layer at low $f_x$, thanks to the additional correction term that models photons with a low number of scattering events/short pathlengths. The modified $\delta \text{-}\mathrm{P}1$ model [Fig. 6(c)] is the most accurate for very thin layers ($d<0.2\mathrm{mm}$), but it seems to over-estimate the contribution of the top layer and becomes less accurate for thickness $d>0.5\mathrm{mm}$. In Fig. 7, we compare the performance of the models by looking at the RMSPE calculated for all thicknesses and three different spectral bands. Because the two scattering agents (${\mathrm{TiO}}_2$ and ${\mathrm{Al}}_2{\mathrm{O}}_3$) have a large difference in the scattering slope, it is possible to choose bands far enough apart, so the ratio ${\mu }_{s_{\text{top}}}^{\prime }/{\mu }_{s_{\mathrm{bot}}}^{\prime }$ is different and we can do an evaluation with different contrasts between the top and bottom layer. In the three bands that we reported (458, 536, and 626 nm), the ratio is equal to approximately 4, 3.5, and 3, respectively.

A first consideration is that all three models have a similar performance for $d>1\mathrm{mm}$, which might be an additional indication that we are exceeding the limits of the two-layers model and measuring only the top layer. Second, both $\delta \text{-}\mathrm{P}1$ models perform better than the SDA in most cases, with the original $\delta \text{-}\mathrm{P}1$ model performing better in the range $0.3\mathrm{mm}<d<0.7\mathrm{mm}$ and our modified $\delta \text{-}\mathrm{P}1$ model performing better in the range $d<0.3\mathrm{mm}$. Third, the RMSPE seems to improve in all models for a smaller ${\mu }_{s_{\text{top}}}^{\prime }/{\mu }_{s_{\mathrm{bot}}}^{\prime }$ ratio (i.e., less contrast between the layers), which is the opposite of what was initially expected. However, when considering that this two-layer approach is based on measurements that assume a homogeneous fluence distribution, it becomes unsurprising that errors in estimating the layer-specific scattering coefficients would increase the more we deviate from this assumption (i.e., for increasing difference in optical properties).

Nevertheless, these preliminary results remain extremely encouraging in the context of the target application (assessment of wound healing), in which we are interested in detecting a difference in ${\mu }_s^{\prime }$ in very thin layers of cells (0.1 to 0.2 mm). In particular, our modified $\delta \text{-}\mathrm{P}1$ model has an overall good performance ($\mathrm{RMSPE}<10\%$), which becomes even better in the physiological range of interest ($\mathrm{RMSPE}<5\%$). The objective of this initial investigation was to evaluate how closely homogenous fluence models could match SFDI measurement of layered media over differing penetration depths.

In Sec. 2.3, we stated our motivations behind the selection of the analytical models of fluence analyzed in this study (i.e., computation speed and ease of implementation). However, any model of light transport suitable for SFDI, which includes spatial frequency information, could theoretically be used without affecting the validity of the two-layer model.^25–27 This allows for future improvement by adopting a model that is capable of representing the fluence of light more accurately in different conditions and for different geometries and optical properties.

Multilayer fluence models in the spatial frequency domain do already exist³⁶ and can be utilized as a second stage optimization, using the results from the homogeneous models as initial estimates of layer thickness and layer specific optical properties.

The current work was aimed at modeling the depth-resolved multi-frequencies SFDI data in the most accurate way possible. However, the method still relies on pre-existent knowledge about the optical and geometrical parameters, which are the ones containing useful information for diagnostic purposes. For the method to be of practical use, an iterative inverse-solving algorithm will be implemented. The algorithm will allow for estimating the scattering parameters and layer thickness from the raw multi-frequencies SFDI measurements, providing the information most useful to aid clinicians in diagnosis. Furthermore, by combining the frequency-dependent depth estimation used in this work with the wavelength-dependent depth estimation from a previous work,²⁰ we will be able to offer a multi-purpose suite of optical tools for the analysis of various skin conditions, ranging from burn wounds to melanoma and everything else in-between.

6.

Conclusions

We have presented an approach to processing SFDI data using different sub-sets of spatial frequencies to obtain datasets that have different depth information. We then made use of this depth-enhanced data to develop and validate a two-layer model of scattering for thin layers, aimed at mimicking the physiology of a healing wound. The model is based on the relative layer contribution to ${\mu }_s^{\prime }$, calculated through an integral function of the light fluence. The performance of three analytical models of fluence was analyzed in the study, but the method itself is model agnostic and can be used with any model of fluence, leaving room for further improvement. The performance of the analytical model themselves looks promising, with the $\delta \text{-}\mathrm{P}1$ models overall working better than the SDA model. The proposed mod-$\delta \text{-}\mathrm{P}1$ model also has an excellent performance for very thin layers, which is especially interesting for the target application of this study, which is to measure layers of skin that are on the order of 0.1 to 0.2 mm.

7.

Appendix A

The initial derivation follows the same procedure seen in the original $\delta$-P1 model.²⁹ By inserting the $\delta \text{-}\mathrm{P}1$ phase function in the Boltzmann transport equation, we obtain the following governing equations:

We can then make the same considerations about the linearity of the medium, giving in response a sinusoidal fluence rate with no phase shift, given as

By introducing Eqs. (14) and (15) in Eq. (11), after the opportune simplifications, we obtain a 1-D differential equation in the depth dimension, given as

From here on, we can continue following the derivations steps found in the appendix of Carp et al.,²⁹ solving the differential equation in a manner similar to the case of planar illumination on a semi-infinite geometry, but substituting the coefficient ${\mu }_{\mathrm{eff}}^2$ with ${\mu }_{\mathrm{eff}}^{\prime 2}$. There are two boundary conditions required. First, because of the conservation of energy, the intensity of the diffuse light field must be zero in regions at a large distance from the source, given as

The second boundary condition is given by the conservation of the diffuse flux component normal to the surface, which is given as

The boundary conditions allow for solving Eq. (16) and finding the diffuse fluence rate, normalized to the incident power $P_0(1-R_s)$, with $R_s$ being the specular reflectance, which is given as

Finally, by adding the collimated (non-diffuse) fluence contribution to Eq. (20), which is given as

8.

Appendix B: Additional Figures

Additional data from the study is reported in this appendix and is shown in Figs. 8Fig. 9–10. In the main text of the paper, the results from the two-layer phantom combinations were presented at 536 nm as a demonstration of the three different models’ performance over increasing spatial frequency sets (Fig. 6). Here, Figs. 8 and 9 also show the same model performance relative to the measured bulk scattering results at 458 and 626 nm, which have a scattering contrast ratio between the layers of approximately 4 and 3, respectively.

Fig. 8

Comparison of the three fluence-based two-layer models (dashed lines) with the data (markers) measured at 458 nm. (a) SDA, (b) $\delta \text{-}\mathrm{P}1$ approximation, and (c) mod-$\delta \text{-}\mathrm{P}1$ approximation.

Fig. 9

Comparison of the three fluence-based two-layer models (dashed lines) with the data (markers) measured at 626 nm. (a) SDA, (b) $\delta \text{-}\mathrm{P}1$ approximation, and (c) mod-$\delta \text{-}\mathrm{P}1$ approximation.

Fig. 10

Absorption coefficient (${\mu }_a$) measured on homogeneous phantoms (solid lines) and two-layered phantoms (markers) at three wavelengths, with the respective scattering contrast of the two-layered phantoms ${\mu }_{s_{\text{top}}}^{\prime }/{\mu }_{s_{\mathrm{bot}}}^{\prime }$ reported in parenthesis: (a) 458 nm (4), (b) 536 nm (3.5), and (c) 626 nm (3).

A particular consideration could be given to the absorption coefficient measurements shown in Fig. 10. Given the manufacturing procedure of the phantoms, we expect a similar absorption in both the ${\mathrm{TiO}}_2$ and ${\mathrm{Al}}_2{\mathrm{O}}_3$ phantoms, which is reflected in the measurements at low $⟨f_x⟩$. However, with increasing $⟨f_x⟩$, we begin to see an increase in ${\mu }_a$ in the multi-layered phantoms, especially the ones with the smallest thickness. We believe that this unexpected deviation is due to a combination of (1) the low sensitivity of absorption at high spatial frequencies (as previously shown in a study from Cuccia¹⁶) and (2) the low signal to noise ratio (SNR) of the spatial frequency dependent reflectance at high spatial frequencies, as can be seen by the large variance of the measurements. This unsettling propagation of error in ${\mu }_a(⟨f_x⟩)$ estimation may be a source of concern regarding the reliability of the proposed models and methods. We found that absorption contributes minimally, if not negligibly, to the fluence estimation at higher spatial frequencies and, subsequently, to this fluence-based approach to interpreting the partial volume contribution of layer-specific scattering. In this higher spatial frequency range, both $⟨f_x⟩$ and scattering become the dominant factors. Nevertheless, further investigation is warranted to determine the exact source of this unexpected variation in ${\mu }_a(⟨f_x⟩)$ as it could benefit future instrumentation design (i.e., impact of a better SNR, dynamic range, and calibration methods in measurements) and model development (i.e., more robust light transport models to reflect increasingly sub-diffusive absorption events).