## 1.

## Introduction

Diffuse optical spectroscopy (DOS) characterizes the optical properties of a medium by measuring the amount of radiative energy remitted from a medium.^{1, 2, 3} Spatial frequency domain imaging (SFDI) is based on DOS principles and uses patterned illumination to determine, with appropriate light transport models, the absorption and reduced scattering coefficients of a turbid medium based on the spatial frequency reflectance function. The method has been described in detail in the literature.
^{4, 5, 6, 7, 8, 9, 10, 11, 12} Instead of measuring the total reflectance,^{1, 2, 3} spatial frequency domain techniques measure the attenuation of specific spatial frequency components of the illuminating pattern as it propagates in a tissue. Illumination by a one-dimensional sine-wave pattern with controlled spatial frequency *f*
_{x} is the simplest way to probe a medium's attenuation of a distinct spatial frequency. The spatial frequency *f*
_{x} can be chosen to have higher sensitivity to a particular physical depth since higher frequencies tend to probe a more superficial volume. Furthermore, higher spatial frequencies have been shown to be more sensitive to the tissue's scattering coefficient. On the other hand, lower frequencies (i.e., *f*
_{x} = 0 mm^{−1}) have been shown to be sensitive to the tissue's absorption coefficient.^{4, 5, 6, 7, 14, 15}

An accurate and efficient model of light transfer through optical media is required to extract quantitative optical properties from SFDI data. Such a model should map the forward and inverse relationships between a set of optical properties and a spatial frequency dependent reflectance function. The radiative transfer equation (RTE) is an accurate description of light propagation in turbid media irradiated with structured light. However, exact solutions to the RTE are known only for a few idealized cases and for planar illumination.^{16, 17} Monte Carlo simulations offer an accurate solution of the RTE and can be adapted to a wide range of multilayered configurations and illumination geometries. However, these are computationally intensive and may not be suitable for real-time applications when immediate estimation of the concentration of tissue constituents such as oxyhemoglobin and deoxyhemoglin concentration are required. The diffusion approximation is frequently used in biomedical optics because it can be a computationally efficient method for estimating light transport in strongly scattering biological tissues. Multiple adaptations of the diffusion approximation exist that account for index mismatch,^{18} multilayered tissue structure, and nondiffuse light sources such as collimated irradiation in plane-parallel media.^{19} This approach has also been used to model light transfer in the spatial frequency domain. However, its applicability is limited to the near-infrared (NIR) since visible light is strongly absorbed by the melanin of the epidermis and the blood in the dermis.^{20, 21} The assumptions of the diffusion approximation become invalid in the visible and UV, with absorption coefficient equaling or exceeding the reduced scattering coefficient.^{22}

This study presents a spatial frequency domain model of two layer media. The forward and inverse model presented here can be used to remove the effects of epidermal absorption from a reflectance signal to facilitate tissue spectroscopy from a large population with varying skin tone. Additionally, it can be used to directly quantify the optical thickness of the epidermis, and thus may be useful in the study of disease mechanisms that thicken or darken the epidermis. To achieve this, an artificial neural network was developed that maps a set of input optical and geometric properties to a spatial frequency dependent reflectance function. An additional artificial neural network was developed that directly solves the inverse problem by mapping a spatial frequency dependent reflectance function to a set of optical and geometric properties, thus avoiding computationally intensive iterative least-squares fitting.

## 2.

## Background

## 2.1.

### Spatial Frequency Domain Reflectance

Spatial frequency domain spectroscopy involves illumination of media with a spatial pattern of the form,

## 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q({x,f_x }) = P_0 ({f}_{x}){\rm cos}(2\pi f_x \,{ x}), \end{equation}\end{document} $$q\left(x,{f}_{x}\right)={P}_{0}\left({f}_{x}\right)\mathrm{cos}\left(2\pi {f}_{x}\phantom{\rule{0.16em}{0ex}}x\right),$$*f*

_{x}is the spatial frequency expressed in mm

^{−1},

*P*

_{0}(

*f*

_{x}) is the incident optical power at the spatial frequency

*f*

_{x}and measured in milliwatts, and

*x*is the spatial coordinate parallel to the mediums surface and measured in millimeters. The remitted energy differs from the illumination pattern due to the optical and geometric characteristics of the sample.

^{4, 5, 6, 7}In fact, the spatial frequency reflectance measured from a turbid medium encodes both depth and optical property information, enabling both quantitation and low resolution depth sectioning of the spatially varying medium optical properties.

^{5}

The depth-sensitivity of spatial frequency domain imaging has been established by several publications.^{5, 8, 9, 10} Cuccia
^{5} demonstrated that spatially modulated illumination facilitates quantitative wide-field optical property mapping and sensitivity to buried heterogeneities in turbid media. They performed spatial frequency domain imaging of tissue simulating phantoms with embedded heterogeneities using 42 spatial frequencies *f*
_{x} between 0 and 0.6 mm^{−1}. They showed that changes in demodulated reflectance at low versus high spatial frequencies are sensitive to the lower versus upper embedded heterogeneities. This work showed that spatial frequency domain spectroscopy can detect contrast between background and heterogeneity, but did not provide a quantitative technique for determining the optical properties of the heterogeneity itself. Konecky
^{8} used spatial frequency domain imaging to detect tube heterogeneities buried in homogeneous tissue simulating phantoms. They measured spatial frequency dependent reflectance from the heterogeneous phantoms at 11 spatial frequencies. They then used an inverse method based on the diffusion approximation to reconstruct tomographic contrast images of the buried tubes. This model, however, may not be applicable if the highly absorbing heterogeneity is close to the surface (i.e., the epidermis of the skin).

## 2.2.

### Artificial Neural Networks

The present work takes a semi-empirical approach toward quantifying heterogeneity in two layered media. Instead of estimating the spatial frequency dependent reflectance with an analytical or approximate model, we performed multiple Monte Carlo simulations and then fit a machine learning algorithm—an artificial neural network—to the output data. An artificial neural network is a data structure that can accurately approximate a nonlinear relationship between a set of input and output parameters from multiple samples of input–output pairs.^{23} Unlike approximate models such as the diffusion approximation, a neural network can be trained to predict tissue reflectance for strongly and weakly absorbing media. Furthermore, a neural network can be trained to directly estimate an inverse relationship between measured tissue reflectance and the tissue's optical properties, thus avoiding iterative techniques such as nonlinear least-squares fitting.

Neural network-based approaches to determining optical scattering and absorption coefficients of biological tissue from tissue and phantom reflectance measurements have been proposed by other investigators.
^{24, 25, 26, 27, 28, 29, 30} For example, Farrell*,* trained a neural network to determine the absorption and reduced scattering coefficient of homogeneous biological tissue from spatially resolved diffuse reflectance measurements at eight source-detector separations.^{27, 28} The authors solved for the spatially resolved diffuse reflectance function for multiple input optical properties using the spatially resolved diffusion approximation. The artificial neural network was then trained to map the spatially resolved reflectance to the set of optical properties. This resulted in a functional inverse relationship between a measurement of reflectance and tissue optical properties without the need to perform iterative least-squares fitting. In fact, the model proposed by Farrell
^{27, 28} was later used by Bruulsema
^{31} to measure changes in the skin's scattering coefficient as a function of changes in blood glucose concentration. Since the underlying function used to generate the training set was the diffusion approximation in a semi-infinite, homogeneous medium, the analysis presented by Farrell is limited to wavelengths in the NIR and for weakly absorbing media that are approximately homogeneous.^{27, 28} On the other hand, Pfefer developed an artificial neural network-based technique for extracting the absorption and reduced scattering coefficients from spatially resolved reflectance measurements of highly absorbing media.^{29} Their approach was similar to that of Farrell,^{27, 28} however, Monte Carlo simulations were used as the underlying photonics model in strongly absorbing media. Additionally, Wang developed a neural network to detect the optical properties of two layer media using spatially resolved reflectance measurements.^{32} They produced a training set from Monte Carlo simulations and then validated the trained neural network on two layer tissue simulating phantoms. They then used their inverse model to detect the absorption and reduced scattering coefficients of two layer media in the ultraviolet and visible ranges, for which absorption coefficient is typically greater than reduced scattering coefficient. However, they reported large prediction errors in all their parameters (absorption and scattering coefficients of the top and bottom layer) that ranged between 0.15 and 1.21 mm^{−1} (20% to 120%).

This study that we have carried out presents a forward and inverse model designed for spatial frequency domain measurements of two layer media. First, the practical limitations of spatial frequency domain imaging are discussed in terms of coupling between top layer absorption coefficient and thickness, and insensitivity of measured reflectance to the top layer's reduced scattering coefficient. Then, a practical description of a two layer tissue model is developed and tested in simulated reflectance spectra from human skin.

## 3.

## Analysis

## 3.1.

### Two Layer Tissue Model

Figure 1 shows the two layer geometry, optical properties, and illumination considered in this study. The superficial layer was illuminated by a collimated light source. The spatial frequency *f*
_{x} of the illumination pattern was considered to range between 0 and 0.25 mm^{−1}. The index of refraction, absorption coefficient, reduced scattering coefficients, and thickness of layer 1 are denoted by *n*
_{1}, μ_{a, 1},
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
, and *d*
_{1}, respectively. The index of refraction *n*
_{1} was assumed to be 1.40 to represent biological tissue.^{33} The absorption coefficient μ_{a, 1} was assumed to range between 0.10 and 2.00 mm^{−1} (Refs.
33, 34, 35, 36, 37, 38, 39). The reduced scattering coefficient
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
was assumed to range between 0.50 and 2.00 mm^{−1}. Finally, the thickness *d*
_{1} was assumed to range between 15 and 150 μm which is typical of many regions of the human body.^{40, 41, 42}

The index of refraction, absorption coefficient, and reduced scattering coefficients of layer 2 are denoted by *n*
_{2}, μ_{a, 2}, and
[TeX:]
${\rm \mu}^ \prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
, respectively. For simplicity, the index of refraction of layer 1 was assumed to be equal to that of layer 2 but not equal to that of air (i.e., *n*
_{1} = *n*
_{2} = 1.40). The absorption coefficient of layer 2, μ_{a, 2}, was assumed to range between 0.01 and 0.20 mm^{−1}, while reduced scattering coefficient
[TeX:]
${\rm \mu}^ \prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
was assumed to range between 0.50 and 2.00 mm^{−1} (Refs.
33, 34, 35, 36, 37, 38, 39).

The spatial frequency dependent reflectance
[TeX:]
$R({n}_1,{d}_1,{\rm \mu }_{{\rm a},1} {\rm,\mu}^ \prime_{{\rm s},1} {\rm,\mu }_{{\rm a},2} {\rm,\mu}^ \prime_{{\rm s},2},{f}_{x})$
$R({n}_{1},{d}_{1},{\mu}_{\mathrm{a},1}{,\mu}_{\mathrm{s},1}^{\prime}{,\mu}_{\mathrm{a},2}{,\mu}_{\mathrm{s},2}^{\prime},{f}_{x})$
was determined with Monte Carlo simulations. Monte Carlo simulation software developed by Wang and Jacques^{43} was used to calculate the radial diffuse reflectance function
[TeX:]
$R({n}_1,{d}_1,{\rm \mu }_{{\rm a},1} {\rm,\mu}^ \prime_{{\rm s},1} {\rm,\mu }_{{\rm a},2} {\rm,\mu}^ \prime_{{\rm s},2},{\rm \rho })$
$R({n}_{1},{d}_{1},{\mu}_{\mathrm{a},1}{,\mu}_{\mathrm{s},1}^{\prime}{,\mu}_{\mathrm{a},2}{,\mu}_{\mathrm{s},2}^{\prime},\rho )$
, where ρ was the radial distance from the simulation's origin. Then, the spatial frequency domain diffuse reflectance function was calculated with the Hankel transform using the method suggested by Cuccia,^{6} namely,

## 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\hskip-5pt R({n}_1,{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu}^\prime_{{\rm s},1},{\rm \mu }_{{\rm a},2},{\rm \mu}^\prime_{{\rm s},2},{f}_{x})\nonumber\\ &&\hskip-5pt\quad = 2\pi \int\nolimits_0^\infty {\rho J_0 ({2\pi f_x })R({n}_1,{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu}^\prime_{{\rm s},1},{\rm \mu }_{{\rm a},2},{\rm \mu}^\prime_{{\rm s},2},{\rm \rho })d\rho },\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{ccc}& & R({n}_{1},{d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{s},1}^{\prime},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime},{f}_{x})\hfill \\ & & \phantom{\rule{1em}{0ex}}=2\pi {\int}_{0}^{\infty}\rho {J}_{0}\left(2\pi {f}_{x}\right)R({n}_{1},{d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{s},1}^{\prime},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime},\rho )d\rho ,\hfill \end{array}$$*J*

_{0}is the 0'th order Bessel function of the first kind. Each simulation was run with 10

^{6}photon packets. For clarity of notation, the spatial frequency domain diffuse reflectance [TeX:] $R({n}_1,{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu}^ \prime_{{\rm s},1} {\rm,\mu }_{{\rm a},2} {\rm,\mu} ^\prime_{{\rm s},2},{f}_{x})$ $R({n}_{1},{d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{s},1}^{\prime}{,\mu}_{\mathrm{a},2}{,\mu}_{\mathrm{s},2}^{\prime},{f}_{x})$ shall be referred to simply as

*R*(

*f*

_{x}).

Figure 2 shows an example of this procedure and illustrates the minimal effect of large changes in
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
on *R*(*f*
_{x}). It shows estimates of the spatial frequency diffuse reflectance *R*(*f*
_{x}) as a function of spatial frequency *f*
_{x} for *n*
_{1} = 1.40, *d*
_{1} between 15 and 150 μm, μ_{a, 1} equaling 1 mm^{−1},
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
between 0.50 and 2.50 mm^{−1}, μ_{a, 2} equaling 0.01 mm^{−1},
[TeX:]
${\rm \mu}^ \prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
equaling 2.00 mm^{−1}, and *f*
_{x} between 0 and 0.25 mm^{−1} predicted by Monte Carlo simulations and Eq. 2. Each bundle of curves was generated by varying
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
between 0.50 and 2.50 mm^{−1} while keeping all other parameters constant. Figure 2 illustrates the weak dependence of *R*(*f*
_{x}) on
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
for the range of *d*
_{1} considered. Large changes in
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
resulted in minimal changes in the reflectance of the two layer medium. Thus, a limitation of the present model is that it is insensitive to
[TeX:]
${\rm \mu}^ \prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
for the range of *d*
_{1} considered. Therefore,
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
was assumed to be equal to 1 mm^{−1} in developing the forward and inverse models presented here.

Additionally, Fig. 3 shows estimates of the spatial frequency dependent reflectance as a function of spatial frequency *f*
_{x} between 0 and 0.25 mm^{−1} predicted by Monte Carlo simulations. It illustrates the strong coupling between μ_{a, 1} and *d*
_{1}. For example, the solid curve indicated by Fig. 3a was generated with μ_{a, 1} equaling 0.1 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, *d*
_{1} equaling 100 μm, μ_{a, 2} equaling 0.01 mm^{−1} and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
equaling 1.7 mm^{−1}. On the other hand, the broken curve that overlays the solid curve was generated with μ_{a, 1} equaling 0.21 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 0.5 mm^{−1}, *d*
_{1} equaling 50 μm, and the same values of μ_{a, 2} and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
. It is apparent that increasing μ_{a, 1} and decreasing *d*
_{1} can produce nearly identical spatial frequency dependent reflectance function. Figure 3 also shows a solid curve indicated by (b) that was generated with μ_{a, 1} equaling 0.2 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 2 mm^{−1}, *d*
_{1} equaling 20 μm, and μ_{a, 2} equaling 0.01 mm^{−1} and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
equaling 1.0 mm^{−1}. Again, a similar spatial frequency dependent reflectance function was generated with μ_{a, 1} equaling 0.1 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, and *d*
_{1} equaling 50 μm for the same values of μ_{a, 2} and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
. Increasing the top layer's thickness has approximately the same effect as increasing its absorption coefficient. This coupling phenomenon was also observed for all
[TeX:]
$\tau^\prime_1 > 0$
${\tau}_{1}^{\prime}>0$
. In fact, for
[TeX:]
$\tau^\prime_1$
${\tau}_{1}^{\prime}$
greater than 4, the top layer was almost completely occluded the bottom layer and the two layer medium effectively became a single layer medium with the optical properties of the top layer. Thus, a limitation of the present model is that the individual values of μ_{a, 1} and *d*
_{1} were not detectable from spatial frequency dependent reflectance. For this reason, this study was limited to determination of the reduced optical thickness
[TeX:]
$\tau^\prime_1$
${\tau}_{1}^{\prime}$
defined as

## 3.2.

### Generating Training Set

The model parameters *d*
_{1}, μ_{a, 1},
[TeX:]
${\rm \mu }_{{\rm a},2},\;{\rm and}\;{\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{a},2},\phantom{\rule{0.28em}{0ex}}\mathrm{and}\phantom{\rule{0.28em}{0ex}}{\mu}_{\mathrm{s},2}^{\prime}$
were sampled according to a uniform distribution,

## 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} f(x) = \left\{ \!\!{\begin{array}{c@{\quad}l} {\displaystyle \frac{1}{{b - a}},} & {a < x < b} \\[8pt] {0,} & {{\rm otherwise}} \\ \end{array}} \right. \end{equation}\end{document} $$f\left(x\right)=\left\{\begin{array}{cc}\hfill {\displaystyle \frac{1}{b-a},}& \hfill a<x<b\\ \hfill 0,& \hfill \mathrm{otherwise}\end{array}\right.$$^{−1}because changing this parameter had little effect of the spatial frequency domain reflectance. Monte Carlo simulations were performed to determine the spatial frequency dependent demodulated reflectance for each sample in the training set.

## 3.3.

### Generating a Validation Set

A validation set was generated to test the performance of the neural network. Ten values for each model parameter *d*
_{1}, μ_{a, 1},
[TeX:]
${\rm \mu }_{{\rm a},2},\;{\rm and}\;{\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{a},2},\phantom{\rule{0.28em}{0ex}}\mathrm{and}\phantom{\rule{0.28em}{0ex}}{\mu}_{\mathrm{s},2}^{\prime}$
were selected along a uniform four dimension grid for a total of 10,000 samples. The number of test samples was chosen such that the entries of the covariance between the model parameters *d*
_{1}, μ_{a, 1}, μ_{a, 2}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
and *R*(*f*
_{x} = 0) for both the training and validation sets were within 1% of each other. This ensured that the test set was statistically close to the training set. As with the training set, the scattering coefficient
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
was equal to 1 mm^{−1}. Then, Monte Carlo simulations were performed to determine the spatial frequency domain diffuse reflectance for each test sample. Training was performed on the training set and assessment of model accuracy was performed with the validation set.

## 3.4.

### Training

A neural network was trained to map a set of input model parameters to a frequency domain diffuse reflectance,

## 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\textit NN_f} ({{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu }_{{\rm a},2},{\rm \mu}^\prime_{{\rm s},2} }) = R(f_x). \end{equation}\end{document} $$\mathit{N}{N}_{f}\left({d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime}\right)=R\left({f}_{x}\right).$$^{−1}, inclusive. The first and second hidden layers had 20 and 5 nodes, respectively. A single layer architecture could not be trained to accurately estimate

*R*(

*f*

_{x}). An architecture with two hidden layers was attempted and the number of hidden nodes in each layer was increased until further addition of nodes stopped improving the performance of the neural network. A hyperbolic tangent sigmoid transfer function was used in the hidden layer and a linear transfer function in the outer layer.

^{23}This transfer function seemed a natural choice due to the appearance of the hyperbolic tangent in many solutions of light transfer problems in multilayered slab geometries.

^{21, 44, 45}The forward model [TeX:] ${\rm NN}_{\rm f} ({{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu }_{{\rm a},2},{\rm \mu} ^\prime_{{\rm s},2} })$ ${\mathrm{NN}}_{\mathrm{f}}\left({d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime}\right)$ was later used as part of a nonlinear iterative least-squares fitting algorithm.

A second neural network was trained to map the spatial frequency domain reflectance *R*(*f*
_{x}) to a set of model parameters
[TeX:]
$\tau _1,\mu _{\rm a,2},\mu ^\prime_{\rm s,2}$
${\tau}_{1},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime}$
, namely,

## 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\textit NN_i} [R(f_x)] = \langle \tau _1,\mu _{\rm a,2},\mu^\prime_{\rm s,2} \rangle, \end{equation}\end{document} $$\mathit{N}{N}_{i}\left[R\left({f}_{x}\right)\right]=\u27e8{\tau}_{1},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime}\u27e9,$$*NN*

_{i}is an inverse relationship. This model was used as a faster alternative to nonlinear iterative least-squares fitting. Four neural networks were trained. Each neural network had 6 input nodes, 2 hidden layers, and 1 output node. The 6 input nodes corresponded to

*R*(

*f*

_{x}) at a discrete spatial frequency between 0 and 0.25 mm

^{−1}. We intentionally chose few spatial frequencies to minimize the time required to acquire spatial frequency domain reflectance data in practice.

^{6}The first and second hidden layers had 20 and 5 nodes, respectively. The output node corresponded to

*d*

_{1}, μ

_{a, 1}, μ

_{a, 2}, or [TeX:] ${\rm \mu} ^\prime_{{\rm s},2}$ ${\mu}_{\mathrm{s},2}^{\prime}$ . A hyperbolic tangent sigmoid transfer function was used in the hidden layer and a linear transfer function in the outer layer.

^{23}

Prior to training, the model parameters and the spatial frequency dependent reflectance were normalized to range between −1 and 1. Training was performed using the MATLAB software package (The MathWorks, Incorporated, Natick, Massachusetts) with Neural Networks and Chemometrics toolbox routines. Each neural network converged in less than 1000 iterations.

## 4.

## Results and Discussion

## 4.1.

### Forward Problem: *NN*
_{f}

Figures 4a, 4b compare estimates for 10,000 validation samples of the diffuse reflectance by the neural network
[TeX:]
$NN_f ({{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu }_{{\rm a},2},{\rm \mu} ^\prime_{{\rm s},2} })$
$N{N}_{f}\left({d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime}\right)$
and Monte Carlo simulation for *f*
_{x} equaling 0 and 0.21 mm^{−1}, respectively, for the validation set with *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1}. The average relative difference between estimates of *R*(*f*
_{x} = 0) by Monte Carlo simulation and the neural network was 0.25%, while the maximum absolute difference was 3.33%. Similarly, the average relative difference between estimates of *R*(*f*
_{x} = 0.21) by Monte Carlo simulation and the neural network was 0.38%, while the maximum absolute difference was 4.34%. Figure 4 illustrates that the neural network generalized the relationship expressed in Eq. 5 from the 50,000 training examples. The goodness of fit between estimates of the reflectance by Monte Carlo simulations and the neural network to a linear model (*y* = *x*) was calculated for 21 values of *f*
_{x} between 0 and 0.25 mm^{−1} and found to be greater than 0.9998 for all cases considered. An r-squared value of unity suggests a perfect linear relationship. It is apparent that the present neural network model exhibits a nearly perfect correlation to Monte Carlo simulations for all *f*
_{x} considered.

## 4.2.

### Inverse Problem with the Least-Squares Fitting

The forward model
[TeX:]
$NN_f ({{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu }_{{\rm a},2},{\rm \mu} ^\prime_{{\rm s},2} })$
$N{N}_{f}\left({d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{a},2},{\mu}_{\mathrm{s},2}^{\prime}\right)$
defined in Eq. 5 was used along with an iterative inverse method to estimate parameters
[TeX:]
${\rm \hat d}_1,{\rm \hat \mu }_{{\rm a},1},{\rm \hat \mu }_{{\rm a},2},{\rm and},{\rm \hat \mu ^\prime }_{{\rm s},2}$
${\widehat{\mathrm{d}}}_{1},{\widehat{\mu}}_{\mathrm{a},1},{\widehat{\mu}}_{\mathrm{a},2},\mathrm{and},{{\widehat{\mu}}^{\prime}}_{\mathrm{s},2}$
from a spatial frequency dependent reflectance *R*
_{ref}(*f*
_{x}). This is done by choosing
[TeX:]
${\rm \hat d}_1,{\rm \hat \mu }_{{\rm a},1},{\rm \hat \mu }_{{\rm a},2},\;{\rm and}\;{\rm \hat \mu ^\prime }_{{\rm s},2}$
${\widehat{\mathrm{d}}}_{1},{\widehat{\mu}}_{\mathrm{a},1},{\widehat{\mu}}_{\mathrm{a},2},\phantom{\rule{0.28em}{0ex}}\mathrm{and}\phantom{\rule{0.28em}{0ex}}{{\widehat{\mu}}^{\prime}}_{\mathrm{s},2}$
to minimize a quadratic cost function,

## 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} L = \sum\limits_{i = 1}^N {[ {NN_f ({\hat d_1,\hat \mu _{\rm a,1},\hat \mu _{\rm a,2},\hat \mu^\prime_{\rm s,2},f_{x,i} }) - R_{{\rm ref}} ({f_{x,i} })} ]^2, } \end{equation}\end{document} $$L=\sum _{i=1}^{N}{\left[N{N}_{f}\left({\widehat{d}}_{1},{\widehat{\mu}}_{\mathrm{a},1},{\widehat{\mu}}_{\mathrm{a},2},{\widehat{\mu}}_{\mathrm{s},2}^{\prime},{f}_{x,i}\right)-{R}_{\mathrm{ref}}\left({f}_{x,i}\right)\right]}^{2},$$*R*

_{ref}(

*f*

_{x, i}) is the spatial frequency dependent reflectance at spatial frequency

*f*

_{x, i}and

*N*is the total number of spatial frequencies. In this study,

*L*was minimized with the Levenberg–Marquardt algorithm implemented in MATLAB. Twenty-one spatial frequencies

*f*

_{x}between 0 and 0.25 mm

^{−1}were used (

*N*= 21). Initial conditions for the iterative minimization were chosen randomly within 50% of the true value of

*d*

_{1}, μ

_{a, 1}, μ

_{a, 2}, and [TeX:] ${\rm \mu} ^\prime_{{\rm s},2}$ ${\mu}_{\mathrm{s},2}^{\prime}$ . Two percent (2%) uniform noise was added to

*R*

_{ref}(

*f*

_{x, i}) to simulate instrumentation error.

The values of *d*
_{1} and μ_{a, 1} could not be found with any certainty with this method. In fact, the converged values
[TeX:]
$\hat d_1$
${\widehat{d}}_{1}$
and
[TeX:]
$\hat \mu _{\rm a,1}$
${\widehat{\mu}}_{\mathrm{a},1}$
were strongly dependent on the initial conditions chosen for minimization while
[TeX:]
$\hat \mu _{\rm a,2}$
${\widehat{\mu}}_{\mathrm{a},2}$
and
[TeX:]
$\hat \mu ^\prime_{\rm s,2}$
${\widehat{\mu}}_{\mathrm{s},2}^{\prime}$
where not. However, the product
[TeX:]
$\hat \tau ^\prime = \hat d_1 (\hat \mu _{\rm a,1} + \hat \mu ^\prime_{\rm s,1})$
${\widehat{\tau}}^{\prime}={\widehat{d}}_{1}({\widehat{\mu}}_{\mathrm{a},1}+{\widehat{\mu}}_{\mathrm{s},1}^{\prime})$
was stable with respect to initial conditions. Figure 5a compares the error between the true value of
[TeX:]
$\hat \tau ^\prime $
${\widehat{\tau}}^{\prime}$
and the value estimated by minimizing *L* in Eq. 7 for *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1}. Figure 5b also shows a histogram of the relative percent error in estimating
[TeX:]
$\hat \tau ^\prime $
${\widehat{\tau}}^{\prime}$
. In this range, the average absolute percent relative error between the input and estimated parameters was 57%. To identify the reason for this high average relative error, we defined a sensitivity of *R*
_{ref}(*f*
_{x}) to τ′ as,

## 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} S(R({f_x }),\tau^\prime) = \frac{{\partial {\textit NN_f} ({{d}_1,{\rm \mu }_{{\rm a},1},{\rm \mu }_{{\rm a},2},{\rm \mu^\prime}_{{\rm s},2} })}}{{\partial \tau ^\prime}}. \end{equation}\end{document} $$S(R\left({f}_{x}\right),{\tau}^{\prime})=\frac{\partial \mathit{N}{N}_{f}\left({d}_{1},{\mu}_{\mathrm{a},1},{\mu}_{\mathrm{a},2},{{\mu}^{\prime}}_{\mathrm{s},2}\right)}{\partial {\tau}^{\prime}}.$$*S*(

*R*(

*f*

_{x}), τ′) as a function of spatial frequency

*f*

_{x}for the case of an optical thick top layer and weakly absorbing bottom layer (τ′ = 0.242, μ

_{a, 1}= 0.05 mm

^{−1}, and μ

_{a, 2}= 1 mm

^{−1}), and optically thin top layer and strongly absorbing bottom layer (τ′ = 0.152, μ

_{a, 1}= 0.20 mm

^{−1}, and μ

_{a, 2}= 1 mm

^{−1}). In both cases, the sensitivity decreases with increasing spatial frequency. While the two curves exhibit similar trends, the sensitivity for an optically thick top layer and weakly absorbing bottom layer is on average 8 times larger than for the opposite case. Consequently, if the top layer thickness

*d*

_{1}was restricted to be greater than 50 μm and μ

_{a, 2}restricted to be smaller than 0.09 mm

^{−1}, the average absolute percent relative error between the input and estimated parameters fell to only 15%.

Figure 7a shows the input and estimated values of μ_{a, 2} for *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1} determined by minimizing *L* in Eq. 7. Figure 8b shows a histogram of the relative percent estimation error in estimating μ_{a, 2}. In this range, the average absolute percent relative error between the input and estimated parameters was 14.7%. Similarly, Fig. 8a shows the input and estimated values of
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
for *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1} determined by minimizing *L* in Eq. 7. Figure 8b shows a histogram of the relative percent estimation error in estimating
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
. In this range, the average absolute percent relative error between the input and estimated parameters was 4.3%.

## 4.3.

### Direct Approach with Neural Network: *NN*
_{i}

Iterative least-squares fitting techniques are effective for determining optical properties from diffuse reflectance measurements. However, the accuracy and computational efficiency of least-squares fitting may be susceptible to initial conditions. A direct mapping between a spatial frequency dependent reflectance and optical properties is desired and was denoted in this study by *NN*
_{i}. A neural network was trained on the same set of data presented in Sec.
4.2. However, six spatial frequencies between 0 and 0.25 mm^{−1} were used as inputs and the parameters
[TeX:]
$\tau ^\prime_1,\mu _{\rm a,2}$
${\tau}_{1}^{\prime},{\mu}_{\mathrm{a},2}$
and
[TeX:]
$\mu ^\prime_{\rm s,2}$
${\mu}_{\mathrm{s},2}^{\prime}$
were used as outputs in an effort to create a direct relationship between measured tissue reflectance and its optical properties

Figure 9a compares the error between the true value of
[TeX:]
$\hat \tau ^\prime $
${\widehat{\tau}}^{\prime}$
and the value estimated with Eq. 6 for *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1}. Figure 9b also shows a histogram of the relative percent error in estimating
[TeX:]
$\hat \tau ^\prime $
${\widehat{\tau}}^{\prime}$
. In this range, the average absolute percent relative error between the input and estimated parameters was 43%. In fact, the direct neural network approach applied to 6 spatial frequencies performed on average 11% better than the least-squares approach applied to 21 spatial frequencies. Additionally, if *d*
_{1} was restricted to be greater than 50 μm and while μ_{a, 2} was kept smaller than 0.09 mm^{−1}, the average absolute percent relative error between the input and estimated parameters fell to only 15% for reasons already discussed.

Figure 10a shows the input and estimated values of μ_{a, 2} for *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1} determined by the direct neural network approach. Figure 10b shows a histogram of the relative percent estimation error in estimating μ_{a, 2}. In this range, the average absolute percent relative error between the input and estimated parameters was 8.0%. Indeed, the mean absolute percent estimation relative error for the direct neural network model is almost half of the error associated with the least-squares fitting method presented in Sec.
4.2.

Figure 11a shows the input and estimated values of
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
for *d*
_{1} between 15 and 150 μm, μ_{a, 1} between 0.10 and 2.00 mm^{−1},
[TeX:]
${\rm \mu} ^\prime_{{\rm s},1}$
${\mu}_{\mathrm{s},1}^{\prime}$
equaling 1 mm^{−1}, μ_{a, 2} between 0.01 and 0.20 mm^{−1}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
between 0.50 and 2.00 mm^{−1} determined by the direct neural network approach. Figure 11b shows a histogram of the relative percent estimation error in estimating μ_{a, 2}. In this range, the average absolute percent relative error between the input and estimated parameters was 3.8%. The neural network performed slightly better in determining
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
than the least-squares fitting method.

Table 1 summarizes the performance of the iterative and direct inverse methods. It shows the mean, mean absolute, and standard deviation of the relative percent difference between input values of
[TeX:]
$\tau^\prime_1$
${\tau}_{1}^{\prime}$
, μ_{a, 2}, and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
and their estimates. The mean error is an indicator of the average bias in prediction of a parameter. For example, the iterative inverse method underpredicts the optical thickness
[TeX:]
$\tau^\prime_1$
${\tau}_{1}^{\prime}$
by 11% while the direct method overpredicts it by 23%. The mean absolute error is an estimate of model accuracy without regard for sign. For example, the direct method exhibits a mean absolute error of 8% in prediction of μ_{a, 2} while the iterative inverse method exhibits a larger error of 14.7%. Finally, the standard deviations reported in Table 1 represent the width of the error distribution. It is apparent that prediction of
[TeX:]
$\tau^\prime_1$
${\tau}_{1}^{\prime}$
exhibits a larger variance than the other parameters and may thus be considered less reliable. Table 1 indicates that the computationally efficient direct method performs as well as the iterative inverse method in predicting μ_{a, 2} and
[TeX:]
${\rm \mu} ^\prime_{{\rm s},2}$
${\mu}_{\mathrm{s},2}^{\prime}$
, but exhibits a larger bias in predicting
[TeX:]
${\rm \tau }_1^{\rm ^\prime }$
${\tau}_{1}^{\prime}$

## Table 1

Mean and standard deviation of the relative percent error between estimates of the model parameters $\tau ^\prime_1,\mu _{\rm a,2}$ τ1′,μa,2 , and $\mu ^\prime_{\rm s,2}$ μs,2′ predicted by minimizing L in Eq. 7 and by the direct inverse method given by Eq. 6 with respect to Monte Carlo simulations. The mean of the absolute relative error is also shown.

Iterative method | Direct method | |||||
---|---|---|---|---|---|---|

Mean | Mean absolute | Std. dev. | Mean | Mean absolute | Std. dev. | |

[TeX:] $\tau ^\prime_1$ ${\tau}_{1}^{\prime}$ | −11% | 57% | 54% | 23% | 43% | 76% |

μ_{a, 2} | 1.6 | 14.7% | 16% | −2.4% | 8.0% | 11% |

[TeX:] ${\rm \mu} ^\prime_{{\rm s},2}$ ${\mu}_{\mathrm{s},2}^{\prime}$ | 0.03% | 4.3% | 4.3% | 0.03% | 3.8 | 4.5% |

## 5.

## Conclusion

This study describes a technique for analyzing spatial frequency dependent reflectance of two layer media. An artificial neural network was used to map input optical properties to a spatial frequency dependent reflectance function of two layer media. Then, iterative fitting was used to determine the optical properties from simulated spatial frequency dependent diffuse reflectance. Additionally, an artificial neural network was trained to directly map spatial frequency dependent diffuse reflectance to sets of optical properties of a two layer media, thus bypassing the need for iteration and significantly reducing the time required for determining tissue optical properties. The present model can be used to determine the optical thickness of a strongly absorbing superficial layer and the absorption and reduced scattering coefficient of a supporting semi-infinite layer. However, the reduced scattering coefficient, absorption coefficient, and thickness of the top layer could not be determined independently.

## Acknowledgments

The authors gratefully acknowledge funding provided by the NIH SBIRs 1R43RR030696-01A1, and 1R43RR025985-01, the NIH NCRR Biomedical Technology Research Center (LAMMP: 5P-41RR01192), the Military Photomedicine Program, AFOSR Grant # FA9550-08-1-0384, and the Beckman Foundation.