2.1.

Spatial Frequency Domain Spectroscopy Instrument and Data Acquisition

Reflectance measurements on a set of layered phantoms spanning a range of optical properties and top layer thicknesses were collected using a SFDS instrument as previously described.¹⁰ Briefly, multispectral spatial frequency-dependent reflectance data was acquired using a broadband light source that is sinusoidally intensity-modulated using a spatial-light modulator and projected onto the turbid phantoms. Fifty-one evenly spaced spatial frequencies were projected onto each sample, ranging from 0 (uniform, planar illumination) to $0.5{\mathrm{mm}}^{-1}$, at $0.01{\mathrm{mm}}^{-1}$ intervals. At each spatial frequency, the illumination pattern was projected at three evenly spaced phase shifts: 0, $2\pi /3$, and $4\pi /3$. This approach allows the use of a simple demodulation scheme to determine the AC magnitude of the spatial frequency-dependent reflectance at a single location.²⁵ The spectral range of the data collected by this SFDS instrument for this particular investigation spans $\lambda =450$ to 1000 nm, at $\sim 1\mathrm{nm}$ spectral resolution (Fig. 1).

Fig. 1

SFDS device and tissue phantom. Broadband optical illumination ($\lambda =450$ to 1000 nm) and sinusoidal spatial modulations with frequencies from 0 to $0.5{\mathrm{mm}}^{-1}$ using phase shifts of 0, $2\pi /3$, and $4\pi /3$ are projected onto a layered siloxane phantom. The resulting images are captured and the strength of the reflected AC modulation is determined.

2.2.

Layered Tissue Phantoms

For this investigation, interchangeable two-layer optical phantom constructs are used to approximate the optical properties of structured tissues such as skin. These two-layer constructs provide a controlled and independently verifiable basis to evaluate the accuracy of our inverse solver. These tissue-simulating phantoms were fabricated using polydimethyl siloxane (PDMS) and follow the general procedure outlined in several previous publications.^11,25,26 In this investigation, the bottom layer phantoms were fabricated to be 3- to 4-cm thick, which may be considered semi-infinite in terms of diffuse reflectance in the visible and near-infrared spectral regions. Specifically, freeze-dried bovine hemoglobin was used to mimic dermal absorption properties across both visible- and near-infrared ranges while titanium dioxide microparticles were used to mimic a range of tissue-relevant scattering properties. Details on the procedure for fabricating these types of hemoglobin-like phantoms are described elsewhere.²⁶ Here, three bottom layer phantoms (labeled phantoms 1, 2, and 3) were fabricated using three distinct concentrations of hemoglobin. Scattering properties are also different between these phantoms but vary by only $\pm 15\%$ [Figs. 2(a)–2(d)].

Fig. 2

Optical properties of the phantoms. (a–d) Depict the reduced scattering coefficient (left axis) and absorption coefficient (right axis) spectra of the top layer (a) and the three possible base layers (b–d) of each phantom. (e) Provides the top layer thickness normalized by the transport mean free paths for both $90\mu \mathrm{m}$ (blue) and $300\mu \mathrm{m}$ (red) thicknesses of the top layer. (f) Depicts the relative scattering to absorption strength of each layer on a log10 scale. The green curve shows the top layer, the black shows base layer 1, the red shows base layer 2, and the blue shows base layer 3.

In this study, we used phantoms of two different top layer thicknesses, 90 and $300\mu \mathrm{m}$, respectively. These were fabricated from the same batch of turbid phantom material, to ensure that the optical properties of each were essentially identical. Naphthol green was selected as the absorbing agent as it provides distinct spectral features from those of the underlying hemoglobin phantoms while also providing absorption across the entire spectral range. Titanium dioxide particles were used as the scattering agent. The resulting absorption and scattering spectra for these two top layer thicknesses are shown in Fig. 2(d). Reference values for the scattering and absorption spectra of the top phantom layer were determined using integrating sphere measurements in conjunction with the inverse adding doubling method,^11,27 whereas bottom layer optical properties were determined using SFDS measurements.¹⁰ Moreover, to understand these layered optical properties vis-a-vis the limitations of the SDA, in Fig. 2(e) we plot the wavelength dependence of the ratio of the reduced scattering coefficient to the absorption coefficient $({\mu }_s^{\prime }/{\mu }_a)$ for each of these layers as well as the top layer thickness normalized to the transport mean free path in Fig. 2(f). Rather than attempt to mimic a specific tissue system, these phantoms were constructed to provide a range of optical parameters and layer thicknesses germane to layered tissue structures (e.g., epithelial tissues) and superficial tissue injury (e.g., burn wounds).²⁸ Figures 2(e)–2(f) show that these properties span many ratios of scattering to absorption (from $\sim 4:1$ to $\sim 100:1$ within a given layer) and layer thickness to $l^{*}$ (from 10:7 to 10:1). Figures 2(a)–2(d) show that either the top or bottom layer can possess the dominant reduced scattering and/or absorption coefficient depending on the wavelength considered.

SFDS measurements were made on six combinations of top and bottom layer phantoms. Each of the three bottom phantoms was measured with either the 90- or $300\mu \mathrm{m}$-thick phantom placed on its top surface. Prior to measurement, all air gaps between the two layers squeezed out through mechanical pressure. SFDS measurements were acquired at 51 spatial frequencies over the range of 0 to $0.5{\mathrm{mm}}^{-1}$ range. Three sets of measurements were acquired from each two-layer phantom configuration. Each set was carried out at a slightly displaced spatial location to average out any potential minor spatial variations in the phantom optical properties. The top and bottom layers were then separated and the top layer was then attached to the another bottom layer phantom. This procedure was performed for all combinations of top and bottom layers, resulting in data from six distinct layered phantom systems. A calibration measurement, using a reference phantom having well characterized optical properties, was performed between each layered phantom^10,24 measurement.

The phantoms were designed with optical property ranges in absorption (${\mu }_a=0.01$ to $0.2{\mathrm{mm}}^{-1}$) and reduced scattering (${\mu }_s^{\prime }=0.5$ to $1.3{\mathrm{mm}}^{-1}$) coefficients that span those of many soft biological tissues in the visible and near-infrared spectral range.²⁸ While the reduced scattering coefficient is dominant over absorption in both layers and across all wavelengths tested, we explore regimes where that dominance is more pronounced in either the top or the bottom layer and systems with top layer thickness far less than the transport free mean path.

2.3.

Spatial Frequency Domain Sensitivities and Design of Multistage Optimization Algorithm

For analysis of this two-layer phantom system, we assume knowledge of layer thickness and attempt to recover four optical properties at each wavelength, namely the absorption and reduced scattering coefficient in each layer. A priori knowledge of layer thickness is a reasonable assumption as existing approaches are available for their estimation using SFDS instrumentation.^24,29 Rather than attempting to fit all parameters simultaneously, we perform a multistage algorithm to isolate and refine the spectra of each optical property. The design of our algorithm is inspired by our prior development of multistaged inversion algorithms^8,9 and the spatial frequency-dependent sensitivity of the measured reflectance to the absorption and reduced scattering coefficients in each layers shown below in Fig. 3.

Fig. 3

Reflectance sensitivities with respect to (a, b) absorption and (c, d) reduced scattering coefficient for (a, c) bottom and (b, d) top layers, respectively, as a function of spatial frequency for a top layer thickness of $300\mu \mathrm{m}$. In each plot, the black line indicates the case where ${\mu }_{s,t}^{\prime }=0.4\times {\mu }_{s,b}$. The red line indicates ${\mu }_{s,t}=0.7\times {\mu }_{s,b}$, the blue indicates ${\mu }_{s,t}={\mu }_{s,b}$, the green indicates ${\mu }_{s,t}=1.3\times {\mu }_{s,b}$, and the magenta indicates ${\mu }_{s,t}=1.6\times {\mu }_{s,b}$.

We present the spatial frequency-dependent reflectance sensitivities with respect to each optical property for a top layer thickness of $300\mu \mathrm{m}$, which is equivalent to the thicker of the layered phantom systems measured in this study. We compute these sensitivities for varying ratios of the reduced scattering coefficient between the top and bottom layers. Specifically, we fix the bottom layer properties at ${\mu }_{s,b}=1{\mathrm{mm}}^{-1}$, ${\mu }_{a,b}=0.1{\mathrm{mm}}^{-1}$ and vary top layer reduced scattering properties while matching the absorption properties with the bottom layer. The variations in top layer scattering that we explore span ${\mu }_{s,t}=[0.4,0.7,1,1.3,1.6]\times {\mu }_{s,b}^{\prime }$ with fixed top layer absorption (${\mu }_{a,t}=0.01{\mathrm{mm}}^{-1}$).

In Figs. 3(a)–3(d), we show the spatial frequency-dependent reflectance sensitivities to (a) bottom layer absorption, (b) top layer absorption, (c) bottom layer reduced scattering, and (d) top layer reduced scattering coefficients. We note that reflectance sensitivity to either top or bottom layer absorption is increrased at the lowest spatial frequencies and decreases by an order of magnitude once spatial frequencies exceed $\sim 0.08{\mathrm{mm}}^{-1}$ for the bottom layer and $\sim 0.13{\mathrm{mm}}^{-1}$ for the top layer. Interestingly, in absolute terms, the reflectance sensitivity to top layer absorption is roughly an order of magnitude smaller than that for bottom layer absorption in this low spatial frequency regime. Moreover, the reflectance sensitivity to absorption in either layer is weakly dependent on the scattering contrast between the layers.

The spatial frequency-dependent reflectance sensitivity to the reduced scattering coefficient in either top or bottom layers shows a distinct peak at a nonzero spatial frequency. The location of the peak sensitivity for the bottom layer is fixed at $0.06{\mathrm{mm}}^{-1}$, whereas the location of this peak for top layer scattering shifts from 0.06 to $0.16{\mathrm{mm}}^{-1}$ as the scattering contrast moves from being top layer dominant (${\mu }_{s,t}=1.6\times {\mu }_{s,b}$) to bottom layer dominant (${\mu }_{s,t}=0.4\times {\mu }_{s,b}$). In absolute terms, the reflectance sensitivity to scattering in the bottom layer can be as much as $5\times$ larger than the top layer sensitivity in this low spatial frequency regime. However, at higher spatial frequencies exceeding $f_x=0.3{\mathrm{mm}}^{-1}$, the reflectance sensitivity to top and bottom layer scattering becomes comparable.

These sensitivity features motivate the design of an inversion algorithm using four stages, each focused on the determination of specific optical parameters. In each stage, we consider measured reflectance values at specific spatial frequencies and 32 equispaced wavelengths spanning $\lambda =450$ to 1000 nm. Using the lsqcurvefit function in MATLAB, we seek to determine optical parameters values that minimize the least squares difference between a ${\mathrm{SHEF}}_N$ computation⁷ for the SFD reflectance and the actual measurements. For the ${\mathrm{SHEF}}_N$ computation, we consider a ninth-order expansion of spherical harmonic functions that provides a good balance between accuracy and computational expense. We will refer to this as a ${\mathrm{SHEF}}_9$ computation for the remainder of the paper. We also chose to consider only 32 discrete wavelengths from the multispectral dataset as it provides a good balance between spectral detail and computational expense. For the ${\mathrm{SHEF}}_9$ computations, we assume a single scattering anisotropy of 0.8 for both layers.²⁸

The details of the four stage process are as follows:

Stage 1: Recovery of layer-specific reduced scattering spectra. In stage 1, we obtain estimates for the reduced scattering coefficient spectra of the top and bottom layers. Because the sensitivity to the top layer relative to the bottom layer increases at larger spatial frequencies [Figs. 3(c) and 3(d)], we choose to analyze separately SFDS data obtained at six spatial frequencies from two different spatial frequency bands ($f_{x1}=[\mathrm{0,0.02},\dots ,0.1]{\mathrm{mm}}^{-1}$ and $f_{x2}=[\mathrm{0.3,0.32},\dots ,0.4]{\mathrm{mm}}^{-1}$). For each frequency band, we analyze the SFDS reflectance spectra using a ${\mathrm{SHEF}}_9$ computation for a homogeneous medium and determine the values for absorption and reduced scattering coefficients that result in a best fit. The presumption here is that the reduced scattering properties obtained from the $f_{x1}$ spatial frequency band will be representative of the bottom layer, whereas those obtained from the $f_{x2}$ spatial frequency band will be representative of the top layer. For simplicity, we assume the spectral dependence of scattering to be governed by an inverse power law²⁸${\mu }_s=A{(\lambda /{\lambda }_0)}^{-b}$, where $A$ is the reduced scattering coefficient for ${\lambda }_0=750\mathrm{nm}$. We discard predictions for optical absorption from this stage.

Stage 2: Recovery of layer specific absorption spectra. In stage 2, we aim to recover estimates for the absorption coefficient spectra of both top and bottom layers. We fix the reduced scattering spectra for the top and bottom layers to those obtained via the Stage 1 analysis of SFDS data in high- and low-spatial frequency bands, respectively. To obtain absorption spectra specific to both top and bottom phantom layers, we consider the SFDS spectral data at spatial frequencies of $f_x=0.01$ and $0.02{\mathrm{mm}}^{-1}$ only. We choose these spatial frequencies as our analysis in Fig. 3 reveals that the reflectance sensitivity to absorption in either layer is maximized at lower spatial frequencies. Moreover, as the reflectance sensitivity to the top layer absorption decays with spatial frequency more gradually than the bottom layer, the SFDS data at $f_x=0.02{\mathrm{mm}}^{-1}$ provide differentially more sensitivity to top layer absorption relative to that obtained at $f_x=0.01{\mathrm{mm}}^{-1}$. Given that the absorption properties of most deep tissue layers are well characterized by a linear combination of the absorption spectra of oxyhemoglobin, reduced hemoglobin, water, and lipid,³⁰ we assume that the spectral shape of the absorption spectra in the bottom layer is known a priori. For the top layer, we adopt a more general approach and do not impose any assumptions regarding the shape of the absorption spectra. We again apply the ${\mathrm{SHEF}}_9$ algorithm, with a two-layer tissue geometry with known top layer thickness (90 or $300\mu \mathrm{m}$). We fit the SFDS data to the ${\mathrm{SHEF}}_9$ predictions by determining the optimal absorption coefficient values for layers 1 and 2 while holding fixed the top and bottom layer reduced scattering coefficient spectra obtained in stage 1.

Stage 3: Refinement of bottom layer reduced scattering spectra. In stage 3, we aim to refine bottom layer reduced scattering spectrum obtained in stage 1. We use SFDS data obtained at $f_x=0.06$ and $0.15{\mathrm{mm}}^{-1}$. These spatial frequencies are chosen as $0.06{\mathrm{mm}}^{-1}$ represents the spatial frequency that is maximally sensitive to bottom layer scattering while $0.15{\mathrm{mm}}^{-1}$ retains significant sensitivity to bottom layer scattering while also effectively eliminating sensitivity to both top and bottom layer absorption. We fix the properties obtained from stage 1 for the reduced scattering spectrum in the top layer and from stage 2 for the top and bottom layer absorption spectra. We drop our constraint for the bottom layer reduced scattering spectra to abide by an inverse power law and perform the fit of the SFDS data to predictions provided by the layered ${\mathrm{SHEF}}_9$ model for each wavelength independently.

Stage 4: Refinement of top layer absorption spectra. The goal of the final stage is to improve the fit of the top layer absorption. As in stage 2, we again use SFDS data from $f_x=0.01$ and $0.02{\mathrm{mm}}^{-1}$ and find the top layer absorption values at each wavelength that produces predictions from our layered ${\mathrm{SHEF}}_9$ model that best match the data. We fix the top layer reduced scattering spectrum to that obtained in stage 1, bottom layer absorption spectrum to that obtained in stage 2, and bottom layer reduced scattering to that obtained in stage 3. For our initial guess, we use the top layer absorption spectrum obtained from stage 2.

3.1.

Stage 1: Recovery of Layer-Specific Reduced Scattering Spectra

Phantom 1: We show the results for phantom 1 in Fig. 4(a) (bottom layer scattering) and 4(d) (top layer scattering). The corresponding errors of these estimates are shown in Figs. 5(a) and 5(d), respectively. For the $300\text{-}\mu \mathrm{m}$-thick top layer, the error in bottom layer reduced scattering is typically $⩽0.1{\mathrm{mm}}^{-1}$ except for $\lambda >450\mathrm{nm}$ where it is $=0.14{\mathrm{mm}}^{-1}$. The error in the recovery of the $300\text{-}\mu \mathrm{m}$-thick top layer reduced scattering coefficient is largest at $\lambda =450\mathrm{nm}$ at $0.256{\mathrm{mm}}^{-1}$ and steadily declines for larger wavelengths. For $\lambda >557\mathrm{nm}$, the error in recovery of the $300\text{-}\mu \mathrm{m}$ thick top layer reduced scattering is $<0.1{\mathrm{mm}}^{-1}$. For the $90\mu \mathrm{m}$ thin top layer phantom, error in bottom layer reduced scattering is worst at $\lambda =628\mathrm{nm}$, where ${\mu }_s$ is underestimated by $0.1{\mathrm{mm}}^{-1}$ with smaller errors occurring at all other wavelengths. The maximal error in recovery of the $90\mu \mathrm{m}$ thin top layer reduced scattering is $0.19{\mathrm{mm}}^{-1}$ at $\lambda =450\mathrm{nm}$ and is reduced at all other wavelengths with errors $<0.1{\mathrm{mm}}^{-1}$ for $\lambda ⩾521\mathrm{nm}$.

Fig. 4

Inversion results. In all graphs, reference data are shown in black, ${\mathrm{SHEF}}_9$ results for a $90\mu \mathrm{m}$ top layer are shown in blue and ${\mathrm{SHEF}}_9$ results for a $300\mu \mathrm{m}$ top layer are shown in red. (a–f) The stage 1 results, in which initial guesses for top and bottom layer scattering spectra are recovered. (a–c) The reference scattering spectra for their respective base layers, and homogeneous ${\mathrm{SHEF}}_9$ results for data taken with the $f_{x1}$ spatial frequency band. (d–f) The reference scattering spectra for the top layer alone and the homogeneous SHEF results using data from the $f_{x2}$ spatial frequency band. (g–i) The top layer absorption spectra fits from stage 2, using results from stage 1 for top and bottom layer scattering. (j–l) The stage 3 results in which recovery of the bottom layer reduced scattering coefficient spectrum is refined. (m–o) The stage 4 results that provide the final result for the top layer absorption spectrum.

Fig. 5

Errors of inversion results. Each subfigure shows the error for thick top layer (red) and thin top layer (blue) phantoms corresponding to the results from the same subfigure shown in Fig. 4.

Phantom 2: The results for phantom 2 are shown in Fig. 4(b) (bottom layer scattering) and 4(e) (top layer scattering.) Error is shown in Fig. 5(b) for top layer scattering and in Fig. 5(e) for bottom layer scattering. For the thick top layer, the error in the bottom layer reduced scattering error is $<0.1{\mathrm{mm}}^{-1}$ for all wavelengths, and $<0.05{\mathrm{mm}}^{-1}$ for all wavelengths apart from $\lambda =450\mathrm{nm}$, where it is $0.076{\mathrm{mm}}^{-1}$. The maximal error in estimation of the $300\text{-}\mu \mathrm{m}$-thick top layer reduced scattering is $0.15{\mathrm{mm}}^{-1}$ at $\lambda =450\mathrm{nm}$, and falls below $0.1{\mathrm{mm}}^{-1}$ for $\lambda ⩾503\mathrm{nm}$ with typical errors of $0.06{\mathrm{mm}}^{-1}$. For the $90\text{-}\mu \mathrm{m}$-thin top layer, a maximal error in the estimation of the bottom layer reduced scattering coefficient is $0.11{\mathrm{mm}}^{-1}$ occurring once again at $\lambda =450\mathrm{nm}$. For $\lambda ⩾503\mathrm{nm}$, this error is $<0.05{\mathrm{mm}}^{-1}$. Top layer reduced scattering error for the $90\text{-}\mu \mathrm{m}$-thin top layer is within $\pm 0.05{\mathrm{mm}}^{-1}$ at all wavelengths.

Phantom 3: Figures 4(c) (bottom layer scattering) and 4(f) (top layer scattering) provide the results from phantom 3. Error is shown in Figs. 5(c) and 5(f), respectively. For the $300\text{-}\mu \mathrm{m}$-thick top layer, error in the recovery of the bottom layer reduced scattering is between $-0.05{\mathrm{mm}}^{-1}$ and 0 at all wavelengths. Top layer reduced scattering error remains within $\pm 0.07{\mathrm{mm}}^{-1}$. For the $90\text{-}\mu \mathrm{m}$-thin top layer, phantom for bottom layer reduced scattering is underestimated by $⩽0.05{\mathrm{mm}}^{-1}$ for all wavelengths. Error in the top layer reduced scattering is within $\pm 0.05{\mathrm{mm}}^{-1}$ for all wavelengths.

3.2.

Stage 2: Recovery of Layer-Specific Absorption Spectra

In stage 2, we determine the magnitude of a predefined spectral shape for bottom layer absorption and an initial estimate for the top layer absorption spectrum. Table 1 shows the results for bottom layer absorption coefficient for each phantom. A recovered coefficient value of $\beta =1$ represents perfect recovery. For all phantoms with the $300\text{-}\mu \mathrm{m}$-thick top layer, error in the recovery of ${\mu }_a$ was $<5\%$. For phantoms 1 and 2, the error is $<1\%$, and for phantom 3, the error is $\sim 4.5\%$. For the $90\text{-}\mu \mathrm{m}$-thin layer phantoms, we recover the $\beta$ value with $<9\%$ error overall. Specifically, the error in the estimation of the bottom layer absorption in phantoms 1, 2, and 3 being 3.6, 6.2, and 8.1%, respectively.

Table 1

Coefficients (β) of the assumed absorption spectrum recovered for each base layer choice and top layer thickness. As the spectrum being fit is the reference data for bottom layer absorption, β=1 indicates perfect recovery.

We show the initial estimates for the top layer absorption in Figs. 4(g)–4(i). It is important to remember that these estimates are provisional; the top layer absorption spectra obtained here will be refined in stage 4. We comment on the accuracy of these estimates at four sets of wavelengths: $\lambda =450\mathrm{nm}$, 450 to 600 nm, 600 to 850 nm, and $>850\mathrm{nm}$. For $\lambda =450\mathrm{nm}$, recovery of ${\mu }_a$ for the $300\text{-}\mu \mathrm{m}$-thick top layer phantom is overestimed between 0.04 and $0.08{\mathrm{mm}}^{-1}$, whereas the $90\text{-}\mu \mathrm{m}$-thin top layer phantoms experience very small errors ($⩽0.01{\mathrm{mm}}^{-1}$). For $\lambda =450$ to 600 nm, the reference data for absorption drops to ${\mu }_a⩽0.07{\mathrm{mm}}^{-1}$, at which point the method experiences difficulties in recovering the correct value, which often results in an extreme underestimation of the absorption. Typically, recovery of ${\mu }_a$ for the $300\text{-}\mu \mathrm{m}$-thick top layer is worse than for the $90\text{-}\mu \mathrm{m}$-thin layer. Here, error for ${\mu }_a$ in $300\text{-}\mu \mathrm{m}$-thick top layer phantoms is between 0.04 and $0.05{\mathrm{mm}}^{-1}$, whereas that of the $90\mu \mathrm{m}$-thin top layer phantoms is near $0.03{\mathrm{mm}}^{-1}$ for phantoms 2 and 3, respectively. Interestingly, recovery of ${\mu }_a$ for the $300\text{-}\mu \mathrm{m}$-thick top layer in Phantom 1 is accomplished with much smaller error than its counterparts with phantoms 2 and 3, whereas the recovery of ${\mu }_a$ in phantom 3 with the $90\text{-}\mu \mathrm{m}$ thin layer has negligible error. For $\lambda =600$ to 850 nm, the reference data for ${\mu }_a$ experience a double hump with values between 0.01 and $0.02{\mathrm{mm}}^{-1}$. Despite some underestimation, this spectral structure is well recovered in the inversion results for each phantom. Errors for ${\mu }_a$ in $300\text{-}\mu \mathrm{m}$ thick top layer phantoms in this region are between 0.02 and $0.05{\mathrm{mm}}^{-1}$, with best performance for phantom 1. Cases with the $90\mu \mathrm{m}$-thin layer have larger error for ${\mu }_a$ at the shorter wavelengths and less at longer wavelengths when compared with $300\text{-}\mu \mathrm{m}$-thick top layer phantoms. These error rates stay in the same general range as those for the $300\text{-}\mu \mathrm{m}$-thick top layer phantoms but tend more toward $0.05{\mathrm{mm}}^{-1}$. For $\lambda >850\mathrm{nm}$, the reference data for ${\mu }_a$ experience a drop off similar to when $\lambda =450$ to 600 nm, falling as $\lambda$ increases until it levels off near $0.05{\mathrm{mm}}^{-1}$. We see a similar underestimation for ${\mu }_a$ in most phantoms in this spectral regime, though here it is more pronounced in the cases with the $90\text{-}\mu \mathrm{m}$-thin top layer. Interestingly, when phantoms 2 and 3 have $300\text{-}\mu \mathrm{m}$ thick top layers, the inversion results are less stable, involving a peak overestimation for ${\mu }_a$ near $\lambda =900\mathrm{nm}$ and a peak underestimation near $\lambda =950\mathrm{nm}$. Error for ${\mu }_a$ here ranges for all phantoms between $\pm 0.05{\mathrm{mm}}^{-1}$. When considering the combined results across all wavelengths, it is important to note that the general shape of the ${\mu }_a$ spectrum is recovered in each phantom, generally within $\pm 0.05{\mathrm{mm}}^{-1}$, and significant underestimations are made when ${\mu }_a<0.075{\mathrm{mm}}^{-1}$. Interestingly, recovery of the absoption spectra for the $90\text{-}\mu \mathrm{m}$-thin top layer phantoms tends to experience smaller absolute errors than for the $300\text{-}\mu \mathrm{m}$-thick top layer phantoms when $\lambda ⩽700\mathrm{nm}$, whereas the opposite is true for $\lambda >700\mathrm{nm}$.

3.3.

Stage 3: Refinement of Bottom Layer Reduced Scattering Spectra

The aim of stage 3 is to improve the initial estimates for bottom layer reduced scattering coefficient spectra obtained in stage 1. We show the results for stage 3 in Figs. 4(j)–4(l), and the relative error rates in Figs. 5(j)–5(l). In all phantoms, error in the recovery bottom layer ${\mu }_s^{\prime }$ was reduced at nearly all wavelengths and with structure of the scattering spectra more faithfully recovered as compared with the results in stage 1. Figure 4(j) shows the results for phantom 1. The reduced scattering coefficient spectrum for the $300\text{-}\mu \mathrm{m}$-thick top layer phantom was recovered with an absolute error of $⩽0.02{\mathrm{mm}}^{-1}$ with the exception of an error of $0.06{\mathrm{mm}}^{-1}$ at $\lambda =450\mathrm{nm}$. The $90\text{-}\mu \mathrm{m}$-thin top layer phantom performed slightly worse, with an underestimation of the reduced scattering coefficient by $0.07{\mathrm{mm}}^{-1}$ at $\lambda =450\mathrm{nm}$ and errors of $\pm 0.04{\mathrm{mm}}^{-1}$ at all other wavelengths. Figure 4(k) shows the results for phantom 2. The absolute error for the reduced scattering coefficient for thick top layer phantom was limited to $\pm 0.05{\mathrm{mm}}^{-1}$ except for $\lambda =450\mathrm{nm}$ where the error was $0.06{\mathrm{mm}}^{-1}$. Recovery of the reduced scattering coefficient for thin top layer phantom was slightly better with error of $\pm 0.05{\mathrm{mm}}^{-1}$. Figure 4(l) shows results in the recovery of the bottom layer reduced scattering for phantoms 3. The performance for both thick and thin top layers is similar with scattering estimated with an error of $⩽0.05{\mathrm{mm}}^{-1}$ across the full spectral region.

3.4.

Stage 4: Refinement of Top Layer Absorption Spectra

The aim of stage 4 is to improve the estimates for top layer absorption obtained in stage 2. We show these results Figs. 4(m)–4(o) and the absolute error of these estimates in Figs. 5(m)–5(o). We consider these results relative to those obtained in ${\mu }_a$ in stage 2.

Phantom 1: The recovered absorption spectra for the top layers in phantom 1 are shown in Fig. 4(m), with error in Fig. 5(m). We obtain improvements in the predicted absorption spectra for both thick and thin layers as compared with the results obtained in stage 2 [shown in Figs. 4(g) and 5(g)]. The overestimation for ${\mu }_a$ found in the case of the thick top layer at $\lambda =450\mathrm{nm}$ is $\sim 0.03{\mathrm{mm}}^{-1}$, whereas the error remains minimal for the case of the thin top layer. We see that error for ${\mu }_a$ is likewise reduced for both the thin and thick top layer cases with improvements between 0 and $0.02{\mathrm{mm}}^{-1}$ in the $\lambda =450$ to 600 nm spectral region. We see greater improvements in ${\mu }_a$ error in the $\lambda =600$ to 850 nm wavelength intervals, particularly for the thick top layer with improved accuracy with ${\mu }_a$ errors generally below $\pm 0.03{\mathrm{mm}}^{-1}$ for the thick top layer and $0.04{\mathrm{mm}}^{-1}$ for the thin top layer. For $\lambda >850\mathrm{nm}$, we obtain slightly worse results until $\lambda =920\mathrm{nm}$ after which the errors obtained in stages 2 and 4 are essentially equivalent. The absorption spectra obtained for phantom 1 in stage 4 improve upon stage 2 for the entire interval of $\lambda =450$ to 920 nm, beyond which the recovery is equivalent. Errors in the recovered top layer absorption are limited to $⩽0.04{\mathrm{mm}}^{-1}$ for $\lambda =450$ to 920 nm.

Phantom 2: Recovery of top layer absorption spectra in phantom 2 is shown in Fig. 4(n), with absolute error shown in Fig. 5(n). For the case with the thick top layer, we see overall improvement in capturing the shape of the absorption spectra with slight reductions in the absolute error. This is seen particularly as a lower (by $0.02{\mathrm{mm}}^{-1}$) overestimation of ${\mu }_a$ when $\lambda =450\mathrm{nm}$ and a less pronounced underestimation in the $\lambda =450$ to 600 nm range (despite having a similar minimum). The double hump of ${\mu }_a$ is once again captured in the $\lambda =600$ to 850 nm spectral range, with smaller (between 0 and $0.01{\mathrm{mm}}^{-1}$ reduction) underestimation for the thick top layer phantom. For $\lambda >850\mathrm{nm}$, both thick and thin top layer phantoms experience a sharp decline as $\lambda$ increases, ending in the near 0 values observed earlier. This is in contrast to the thick top layer performance from stage 2, which involved an overestimation of ${\mu }_a$ when $\lambda =880$ to 920 nm. Errors remain similar in absolute value to those of stage 2 but are more uniform across wavelength and better represents the overall shape of the reference absorption spectrum.

Phantom 3: Recovery of top layer absorption spectra in phantom 3 is shown in Fig. 4(o), with absolute error shown in Fig. 5(o). For the thick top layer phantom, the initial overestimation at $\lambda =450\mathrm{nm}$ is actually worse, increasing to $0.05{\mathrm{mm}}^{-1}$ as compared with $0.04{\mathrm{mm}}^{-1}$ in stage 2. While the absolute error is worse in many spectral regions, the recovered spectral shape is better captured as compared with stage 2, and we underestimate the absorption by only $0.05{\mathrm{mm}}^{-1}$ in the case of the thick top layer phantom at $\lambda =588\mathrm{nm}$, with a smaller underestimation ($0.03{\mathrm{mm}}^{-1}$) for the thin top layer phantom. Absolute error in absorption stabilizes to within $\pm 0.04{\mathrm{mm}}^{-1}$ for $\lambda =605$ to 795 nm. For longer wavelengths, we recover a more accurate spectral shape for the case of the thick top layer and comparable results for the case of the thin layer as compared with stage 2. This manifests in the same way it does for phantom 2, with an elimination of the overestimation found in stage 2, once again resulting in a more accurate recovery of the spectral structure of ${\mu }_a$.