1.

Introduction

Quantitative spatial frequency domain imaging (SFDI) for separation of the absorption (${\mu }_{\mathrm{a}}$) and the reduced scattering coefficient (${\mu }_{\mathrm{s}}^{\prime }$) gained importance in biomedical imaging over the last years. It facilitates the investigation of different kinds of tissue changes and biological processes, e.g., blood perfusion, neuronal death, and brain cellular composition changes.^1–3 To eliminate specular and surface reflections in SFDI, crossed linear polarizing filters are generally used.⁴ However, as shown by Wiest et al.,⁵ the use of crossed polarizers causes significant errors in ${\mu }_{\mathrm{s}}^{\prime }$ and ${\mu }_{\mathrm{a}}$ when using a nonpolarized model for solving the inverse problem. As it is shown in this work, surface scattering $R_S$ leads to an offset in the spatial frequency resolved reflectance $R_{\mathrm{SFD}}(F)$. This offset mostly affects high spatial frequencies $F$, as for those $R_{\mathrm{SFD}}(F)$ is steadily decreasing^6,7 and therefore, the ratio $\frac{R_s}{R_{\mathrm{SFD}}(F)}$ is increasing. Correcting for this offset without any more information than $R_{\mathrm{SFD}}(F)$ is hardly possible, because the subdiffusive scattering properties (e.g., phase function) and the surface scattering affect the spatial-frequency resolved reflectance in a similar way. Bassi et al.⁸ were the first to investigate the phase shift due to turbid media, but the phase information was not used for solving the inverse problem. By now, the phase information was only used for determining the topography of the sample in order to correct for changes of the reference intensity.⁹

First, results associated to this work were already presented on the European Conference on Biomedical Optics in 2017.¹⁰ However, the modified solution of the radiative transfer equation (RTE) is now presented in more detail and new results of measurements are shown. The modified solution makes it possible to correct for surface scattering by solving the inverse problem not only taking the reflectance $R_{\mathrm{SFD}}(F)$ but also the spatial shift or phase ${\varphi }_{\mathrm{SFD}}(F)$ into account. It is shown that this phase is significantly influenced by surface scattering and thus appropriate for separating volume and surface scattering. The theoretical results were confirmed by SFDI measurements using specially prepared solid phantoms.

2.

Theory

In order to model the light propagation in turbid media, an analytical solution of the RTE for semi-infinite media is applied within this work. The analytical solution of the RTE is based on a series expansion in rotated spherical harmonics of order $M$. A detailed derivation can be found in literature.^11,12 The solution of this derivation gives the complex reflectance:

The reflectance can be illustrated [compare Fig. 1(a)] as the amplitude of a sinusoidal illuminated pattern, which is being attenuated subject to the optical properties of the sample. For measurement of the reflectance, the sample is illuminated with a periodical intensity pattern, e.g., a collimated sinusoidal intensity of the kind:

Fig. 1

(a) Schematic illustration of the radiance inside a turbid medium, which is obliquely illuminated at ${\theta }_0$ by a sinusoidal intensity pattern. (b) Schematic illustration of the light being detected in a SFDI setup. A certain part $R_V$ of the detected light is coming from the volume (chessboard) and a second part $R_S$ is due to surface scattering (black solid), which both are the normalized integrals of the radiance $L_V$ and $L_S$, respectively.

Experimentally, it is impossible to perfectly align the sample. For this reason, the distance $z$ of misalignment has to be included into the model in order to correctly fit the phase. This is why the spatial shift illustrated in Fig. 1(a) in total reads

If not stated otherwise, whenever the phase is fitted the misalignment, $z$ is always used as additional fit parameter. But as this is not crucial for separation of volume and surface scattering, the misalignment will not be discussed any further within this work.

For samples with an optically smooth surface, which has a $\text{Scratch}\text{-}\mathrm{Dig}<80$,¹³ the total reflection ${\tilde{R}}_{\mathrm{S}}$ is given by Fresnel’s equation according to

This solution describes the light propagation in a semi-infinite turbid medium exact as long as the direct reflex does not fall into the aperture of the detection system.

For a large number of practically relevant cases, the requirement of an optically smooth surface is not fulfilled. In the case of rough surfaces, the detected light is always a superposition of light reflected from the turbid volume and light that is diffusely reflected at the surface. In order to describe the ratio of light, which is scattered at the surface into a certain direction $\overrightarrow{s}({\theta }_S,{\varphi }_S)$ with spherical coordinates ${\theta }_S$ and ${\varphi }_S$, the so-called bidirectional reflectance distribution function (BRDF):

A simple assumption that can be made for the angular distribution of the light scattered at the surface is a constant BRDF. However, as Kienle and Foschum¹⁴ showed such Lambert surfaces do not exist in practice, thus, more precise models have to be applied. For instance, Sun¹⁵ utilize a statistical model of the surface for approximating the angular surface scattering function.

However, in order to only correct for surface roughness when solving the inverse problem, the exact knowledge of the BRDF is even not necessary, as the detected intensity

Within this work, the Lambert approximation ${\rho }_{\mathrm{BRDF},S}=\text{constant}$ was used, for which Eq. (17) together with Eq. (19) yields

It is moreover assumed that the surface scattering can be treated as a small perturbation to the Fresnel equations, which means that the light transmitted both into and out of the sample is disturbed much less by the rough surface than by volume scattering. This assumption also implies that the surface BRDF has no spatial dependency and surface scattering can be treated as local problem. In physical manner, this assumption is fulfilled, if both, the RMS surface roughness $R_q$ and correlation length $l_c$ are much smaller than the mean free path $1/{\mu }_s$ of the volume scattering. It should also be mentioned that the model is restricted to statistical surfaces for which the Kirchhoff approximation (KA) is valid. According to Ulaby et al.,¹⁶ when considering a stationary isotropic Gaussian surface, the KA can be applied, if

In the case of a sinusoidal irradiation [compare Fig. 1(a)], the detected intensity $I_{q_i}=I_{V,q_i}+I_{S,q_i}$ can be understood as a superposition of surface and volume scattered light. For clarification, the difference between these two parts is schematically illustrated in Fig. 1(b), which shows the part of light that was diffusely reflected at the surface [black solid part in Fig. 1(b)] and light that entered the turbid medium and was then scattered back into the objective lens [chessboard in Fig. 1(b)]. The sum

The amplitude measured for an angular spatial frequency $q_i$, thus, yields $\frac{I_0}{2}{R}_{\mathrm{SFD}}(q=q_i)$, which in future is denoted as AC reflectance. The part $\frac{I_0}{2}(R_S+R_V(q=0{\mathrm{mm}}^{-1}))$ corresponds to the intensity, which would be measured utilizing an homogeneous illumination of intensity $I_0/2$, this is furthermore denoted as DC part. For modulated illuminations with $q\ne 0{\mathrm{mm}}^{-1}$, this DC part can be seen as the offset of the sine along the ordinate. According to Fig. 1(a), the total shift along $x$ becomes $x_{\varphi }={\varphi }_{\mathrm{SFD}}(q=q_i,z)/q_i$, where ${\varphi }_{\mathrm{SFD}}$ corresponds to the total measured phase shift, originating both from the turbidity of the medium and the surface reflection. As an example, Fig. 2 shows a calculation of the SFD reflectance $R_{\mathrm{SFD}}$ and phase ${\varphi }_{\mathrm{SFD}}$ versus spatial frequency $F=q/2\pi$. The optical properties for this calculation were ${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}=4.00{\mathrm{mm}}^{-1}$, and $n_{\mathrm{med}}=1.4$, whereby a Henyey–Greenstein phase function with $g=0.75$ was utilized resulting in ${\mu }_{\mathrm{s}}^{\prime }=1.00{\mathrm{mm}}^{-1}$. The total reflectance and phase are shown as green solid lines in Fig. 2 and those are composed of the part coming from the volume $R_V$ (blue dashed) and the surface $R_S$ (orange dotted). For this calculation, the surface scattering factor $R_S$ was exemplarily chosen to be 0.05, as for skin it is expected to obtain comparable values. The influence of surface scattering regarding $R_{\mathrm{SFD}}$ and the changes of the phase can clearly be seen.

Fig. 2

Theoretical calculations of (a) the reflectance $R_{SFD}$ and (b) the phase ${\varphi }_{SFD}$ for different spatial frequencies $F=q/2\pi$ and an irradiation angle of ${\theta }_0=45°$ are shown. The optical properties are ${\mu }_a=0.01{\mathrm{mm}}^{-1}$, ${\mu }_s=4.00{\mathrm{mm}}^{-1}$, and $n_{\mathrm{med}}=1.4$, whereby a Henyey–Greenstein phase function with $g=0.75$ was utilized resulting in ${\mu }_s^{\prime }=1.00{\mathrm{mm}}^{-1}$. The blue dashed curve denoted with $R_V$ or rather ${\varphi }_V$ shows the part which comes from volume scattering. The dashed blue line named $R_S$ or rather ${\varphi }_S$ represents the surface scattered part, for which the surface scattering correction parameter $R_S=0.05$ was exemplarily chosen. The solid green line represents the totally detected reflectance $R_{SFD}$ composed of the two parts $R_V$ and $R_S$.

In order to be able to neglect the influence of the phase function ${\rho }_V(\theta )$ on the determination of the other optical properties, the exact phase function of the used scatterers has to be taken into account. As the solution of the RTE is a series expansion in spherical harmonics, whereby in the following, $M$ denotes the order of this expansion, the phase function likewise has to be expanded in terms of Legendre polynomials. This is why the solution of the transport equation takes $i\in [0;M+1]$ expansion coefficients:

3.

Materials and Methods

All SFDI measurements within this work were carried out with the spatial frequency domain setup described by Bodenschatz et al.¹⁷ It consists of a tungsten polychromatic light source, a filter wheel equipped with eight bandpass color filters, a digital micromirror device (0.7 XGA VIS, Discovery 4100 Development Kit, Vialux, Germany), and a cooled charge coupled device (CCD) camera (QSI640, USA). The sinusoidal intensity patterns with possible spatial frequencies of $0{\mathrm{mm}}^{-1}\le F\le 1{\mathrm{mm}}^{-1}$ are projected under an oblique angle of ${\theta }_0=35\mathrm{deg}$ onto the sample and diffuse reflected light from a $38{\mathrm{mm}}^2\times 38{\mathrm{mm}}^2$ region is imaged by an objective lens ($\mathrm{NA}\approx 0.1$) onto the CCD chip.

By projecting three sine patterns of $J$ and different spatial frequencies $F_i$, each phase shifted by 0, $2\pi /3$, and $4\pi /3$, images with intensity values $I_{n,(u,v)}$ $(F_i)$ for $n\in [\mathrm{1,2},3]$ at pixel $(u,v)$ are acquired. Demodulating each pixel^18,19 yields the DC intensity as follows:

Inhomogeneities in the intensity of the incident light field as well as the point spread function of the illumination system cause artifacts in the demodulated signal, this is why additionally, a reference phantom with known reflectance $R_{\mathrm{SFD},\mathrm{ref}}(F_i)$ and phase ${\varphi }_{\mathrm{SFD},\mathrm{ref}}(F_i)$ has to be measured. With help of this reference measurement, the absolute SFD reflectance is

Note that hereinafter, omitting the declaration $(u,v)$ indicates an averaging over a certain region of the image for which the error is given by the standard deviation of the mean value. All measurements were carried out with the SFDI setup described above at spatial frequencies $F_1=0.043{\mathrm{mm}}^{-1}$, $F_2=0.086{\mathrm{mm}}^{-1}$, $F_3=0.130{\mathrm{mm}}^{-1}$, $F_4=0.173{\mathrm{mm}}^{-1}$, $F_5=0.216{\mathrm{mm}}^{-1}$, and $F_6=0.432{\mathrm{mm}}^{-1}$. Afterward, the model introduced in Sec. 2 was fitted to the data using a nonlinear least squares DOGLEG²⁰ algorithm, which was taken from the library CERES SOLVER.²¹ An example of such a fit is given in Sec. 4.1.

The utilized turbid samples (shown in Fig. 3) were all made from epoxy resin (Crystal Resin Giessharz glasklar, SuK Hock GmbH, Germany) with titanium dioxide added as scatterers. A precise description of the production process and a detailed characterization (e.g., absorption and scattering coefficient) of these kinds of solid-state phantoms was given by Krauter et al.²² Nine different phantoms were produced, from which three are nonturbid with ${\mu }_s^{\prime }\approx 0{\mathrm{mm}}^{-1}$ [Fig. 3(a), pa0], the next three have reduced scattering coefficients of ${\mu }_s^{\prime }\approx 1{\mathrm{mm}}^{-1}$ [Fig. 3(a), pa1] and another three have ${\mu }_s^{\prime }\approx 4{\mathrm{mm}}^{-1}$ [Fig. 3(a), pa4]. Within one group of phantoms with same reduced scattering (e.g., pa4), each of the phantoms were grinded with sand of different grain sizes resulting in three different surface roughnesses [Fig. 3(a), s1 to s3]. These surfaces were then analyzed with a contact profilometer (LV50E, Hommelwerke, Germany), which measures a topography profile $y(x)$ along a line of $l_x=4.8\mathrm{mm}$ length by scanning a diamond stylus (tip radius $5\mu \mathrm{m}$, tip angle 90 deg) in contact across the surface. Three of such profiles were measured at different positions for each of the phantoms pa0-s1, pa0-s2, and pa0-s3. The autocorrelations of the measured surface profiles revealed a Gaussian distributed profile for each of the surfaces; hence, the correlation length $l_c$ could be determined by fitting the normalized autocorrelation with the correlation model:

Fig. 3

(a) Nine epoxy resin phantoms with different volume and surface scattering are shown. From top to bottom (pa0, pa1, pa4), the volume scattering increases from ${\mu }_s^{\prime }\approx 0\mathrm{mm}$ to ${\mu }_s^{\prime }\approx 4\mathrm{mm}$ and from left to right (s0, s1, s2), the surface roughness of the phantoms increases, as indicated by the surface scattering parameter $R_s$, which has been introduced in Sec. 2. This increase is clearly visible for the first row, for which the different shades of gray are only due to the various surface roughness. (b) A epoxy resin phantom, which contains a rough structure on its surface, is shown. The structure is almost invisible for the naked eye, this is why for convenience, the structure is additionally outlined in blue. The red box displays the region, for which the SFDI was done and the investigated model was fitted.

Table 1

The table gives the root mean square parameter Rq and maximum height of profile RZ calculated according to the ISO definition for the three different surfaces s1, s2, and s3 the measured surface profile. The table also shows the correlation length lc, which could be determined form the autocorrelation of the profile, when considering a stationary isotropic Gaussian surface. The last column of the table shows that the KA is, according to Ulaby et al.,16 valid for all of the investigated surfaces.

Thus, a three times, three varieties of phantoms with different surface roughness (s1 to s3) and reduced scattering coefficient (pa0, pa1, and pa4) were investigated. No additional absorbers were added to the phantoms, hence, the absorption is only that of the pure epoxy resin $[{\mu }_{\mathrm{a}}(600\mathrm{nm})\approx {10}^{-3}{\mathrm{mm}}^{-1}]$, for which the absorption spectrum in detail is given by Krauter et al.²²

Additionally, one phantom with homogeneous optical properties of ${\mu }_{\mathrm{s}}^{\prime }\approx 1.05{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{a}}\approx 0.44{\mathrm{mm}}^{-1}$ for the wavelength $\lambda =700\mathrm{nm}$ was produced. The surface was then treated by sand blasting together with a template such that a certain region of the surface got rough and the remaining surface stayed smooth. Both the smooth and the rough surface were then characterized with the same contact profilometer mentioned before. This profile measurement revealed $R_q=(0.08\pm 0.02)\mu \mathrm{m}$ and $R_Z=(1.4\pm 0.3)\mu \mathrm{m}$ for the smooth region A, outlined in Fig. 6(a). The measurement of the rough region B, again displayed in Fig. 6(a), yielded $R_q=(1.8\pm 0.3)\mu \mathrm{m}$ and $R_Z=(11.6\pm 0.5)\mu \mathrm{m}$. A comparison of these values with the parameters given in Table 1 shows that the smooth and the rough surfaces are comparable with the surfaces s1 and s2, presented before. An image of this phantom is shown in Fig. 3(b) with blue lines surrounding the regions of the rough surface, which due to the volume scattering can hardly be seen by the naked eye. The red outlined box in Fig. 3(b) displays the region from which SFDI data were obtained and later investigated by means of the optical properties and surface roughness.

4.

Results and Discussion

Two different sorts of phantoms were investigated within this work. In the following, the measurement results for these phantoms are presented and discussed.

4.1.

Homogeneous Surface Roughness Phantoms

In this section, the results of the measurements with the phantoms shown in Fig. 3(a) are presented. SFD reflectance measurements of these phantoms, as described in Sec. 3, were carried out in order to obtain an averaged $R_{\mathrm{SFD}}$ and ${\varphi }_{\mathrm{SFD}}$. The top plots in Figs. 4(a) and 4(b) exemplarily show the reflectance and phase for the phantom pa1-s2. The green markers indicate the measured data, which was fitted by two models, one which does not take surface scattering into account (blue solid line) and the new model that considers the surface scattering parameter $R_S$ (orange dashed line). The bottom plots in Figs. 4(a) and 4(b) show the relative deviation $(\text{measurement}-\text{fit})/\text{measurement}$ between measurement and fitted model respectively. It can be seen that the modified model fits the data far better than the model that does not contain the surface scattering parameter. The deviations for the model without $R_S$ are almost 1 for the reflectance and 20 for the phase. In comparison, the fit of the modified model, which takes $R_S$ into account, has a maximum deviation of 0.10 for the reflectance and 0.65 for the phase. This is more than a magnitude smaller for both, reflectance and phase, compared to a fit without surface scattering parameter $R_S$.

Fig. 4

The green markers indicate the measured (a) SFD reflectance and (b) phase of the phantom pa1-s2 [see Fig. 3(a)], which were fitted by the model introduced in Sec. 2. The solid blue line shows the results for a fit without taking the surface scattering parameter $R_S$ into account, whereas the dashed orange line displays the results with $R_S$ being considered. The relative deviation between measured data and fitted model is given in the subplot at the bottom. The plots in panel (c) show the predicted against the known (from pa0) surface scattering parameter $R_S$ for four different wavelengths. The circular markers represent the fitted $R_S$ parameters for the group of low scattering phantoms pa1 (${\mu }_s^{\prime }\approx 1{\mathrm{mm}}^{-1}$) and the quadratic markers the fitted surface scattering values for the group pa4 of phantoms with a reduced scattering coefficient of about $4{\mathrm{mm}}^{-1}$. The dashed gray line indicates the target values with a $\pm 3\%$ trust region, highlighted in brighter gray.

Fitting the modified theory yields besides the optical properties ${\mu }_{\mathrm{s}}^{\prime }$ and ${\mu }_{\mathrm{a}}$, the new parameter $R_S$, which as discussed in Sec. 2, is a measure of the surface scattering. In order to validate the correctness of $R_S$, SFD images of the three phantoms pa0-s1, pa0-s2, and pa0-s3 [compare Fig. 3(a), top row] were taken from which $R_{\mathrm{SFD}}$ could be demodulated according to Eq. (29) and as $R_V=0$ for these phantoms relation, Eq. (25) yields

Eq. (33)

$$R_S=R_{\mathrm{SFD}}.$$

These $R_S$ values were taken as reference and are denoted $R_{S,\mathrm{ref}}$. The trust region of $\pm 3\%$ for these values was estimated by multiple measurements of phantoms form group pa0. A validation of the correctness of $R_S$ after solving the inverse problem for the phantom groups pa1 and pa4 [compare Fig. 3(a) middle and bottom row], which both contain volume scattering, is given in Fig. 4(c). The green circles display the result for the phantom group pa1 against the reference values $R_{S,\mathrm{ref}}$ and the green squares indicate the results for the phantom group pa4, respectively. The dashed gray line indicates the target values with $\pm 3\%$ trust region, which is highlighted in brighter gray. The calculated $R_S$ values for both groups pa1 and pa4 are in average about 1% higher than expected, which is probably because of small variations in the measured reflectance values of the highest spatial frequency $F_6=0.432{\mathrm{mm}}^{-1}$. Theoretical investigations showed that especially high spatial frequencies are important for the correct separation of volume and surface scattering.

Solving the inverse problem without surface scattering implies a significant overestimation of the reduced scattering and absorption coefficient. This behavior is shown in Figs. 5(a) and 5(b) for each of the two phantom groups pa1 with ${\mu }_{\mathrm{s}}^{\prime }\approx 1{\mathrm{mm}}^{-1}$ and pa4 with ${\mu }_{\mathrm{s}}^{\prime }\approx 4{\mathrm{mm}}^{-1}$. The plots in Fig. 5 display ${\mu }_{\mathrm{s}}^{\prime }$ against the different surface roughnesses $s1$, $s2$, and $s3$ quantized by the known surface scattering parameter $R_S$. The dashed gray lines indicate the known reduced scattering coefficients with a trust region of $\pm 5\%$ for different wavelengths. The trust region of ${\mu }_{\mathrm{s}}^{\prime }$ was estimated by the production accuracy of the epoxy resin phantoms according to Krauter et al.²². The green circles and squares in Figs. 5(a) and 5(b) are the results obtained with the modified theory, which takes surface scattering into account. They all agree within $\pm 0.25{\mathrm{mm}}^{-1}$ with the reference, which is given as gray dashed line. The red dots and red stars in Figs. 5(a) and 5(b) indicate the values calculated without taking surface scattering into account. These values show a systematic increase of ${\mu }_{\mathrm{s}}^{\prime }$ subject to surface scattering of the sample. The investigations show that without considering surface scattering in the utilized theory, the volume scattering is systematically overestimated.

Fig. 5

The plots in (a) and (b) show the fitted reduced scattering coefficients ${\mu }_{\mathrm{s}}^{\prime }$ against the known surface scattering parameter $R_S$ evaluated with different models, whereas in panel (a), the results for the group pa1 of weakly volume scattering phantoms are displayed, and in panel (b), the results of the highly scattering phantoms pa4 are shown. Green circles and squares have been evaluated with the new model, which contains $R_S$. The red dots and red stars as well as the black crosses have been determined with a model, which does not take surface scattering into account. Note that the SFDI data for the black crosses were acquired with crossed linear polarizing filters between irradiation and detection. The dashed gray line indicates the known values with a trust region of $\pm 5\%$ relative deviation.

In literature, e.g.,^9,24,25 it is often stated that the use of crossed linear polarizers corrects for diffusely surface scattered light. Thus, in addition to the measurements with unpolarized light measurements with crossed linear polarizers between irradiation and detection were carried out both for the reference and the actual sample. The SFDI data obtained from these measurements were then evaluated with the solution of the unpolarized RTE. The black crosses in Figs. 5(a) and 5(b) mark ${\mu }_{\mathrm{s}}^{\prime }$ found by solving the inverse problem. It appears that the crossed polarizers only reduce the influence of the surface scattered light a bit, but the overestimation of ${\mu }_{\mathrm{s}}^{\prime }$ due to a rough surface is still present. Another method for reducing the surface scattering is to emerge the samples in a liquid (e.g., water), but as no analytical solution for a nonturbid layer on top of the actual scattering sample was available and moreover the surface scattering would still be present in a reduced kind, emerging the sample in a liquid was not further investigated within this work.

4.2.

Structured Surface Roughness Phantom

In this part, the modified theory is applied for SFDI measurements of the spatially structured phantom shown in Fig. 3(b), which has both known scattering and absorption properties. Moreover, as already mentioned in Sec. 3, the glossy surface of the phantom has a structure with rough surface in the middle. After SFDI data were acquired as described in Sec. 3, the inverse problem was solved for both the modified theory which contains $R_S$ and a model, which does not correct for surface scattering. The calculated $R_S$ map for a wavelength of 700 nm can be seen in Fig. 6(a), for which the structured rough surface is clearly contrasted compared to the glossy surface surrounding the structure. A box plot of $R_S$ for the two outlined regions (A) and (B), glossy and rough, respectively, is given in Fig. 6(b). The values of the glossy surface A lie at around $R_S\approx 0.8\%$ and for the rough surface at about $R_S\approx 4.6\%$ and are, therefore, in the expected range. For the glossy surface, a roughness parameter of zero was striven for, but due to the hardening process of the initially liquid epoxy resin, an optical exactly smooth surface cannot be achieved. For this reason, the obtained surface scattering parameter $R_S$ for the glossy surface of the structured phantom seems reasonable. The rough surface was obtained by sand blasting with a grain size comparable to that of the particles utilized for grinding the phantoms with surface s2, which were discussed in the previous Sec. 4.1 and are shown in Fig. 3(a). The investigations discussed in the previous Sec. 4.1 yield $R_S\approx 5.5\%$ for such surfaces, which also agrees to the $R_S$ values obtained for the rough sandblasted structure.

Fig. 6

(a) The fitted $R_S$ parameter inside the red outlined region given in Fig. 3(b). is shown. (b) The box plot shows the averaged $R_S$ value of the outlined regions A and B.

Maps of the reduced scattering coefficient of the structured phantom are shown in Figs. 7(a) and 7(b), respectively. The inverse problem for the map Fig. 7(a) was calculated under consideration of $R_S$, whereas map Fig. 7(b) was determined without taking surface scattering into account. As the volume scattering and the reduced scattering coefficient of both regions should be the same, no structure should be present at all in Fig. 7(a). However, in both maps, the structure can clearly be seen, but compared to the map in Fig. 7(b), the reduced scattering coefficient in Fig. 7(a) is only different at the boundaries between smooth and rough surface. The evaluation of the SFDI data without $R_S$ yields ${\mu }_{\mathrm{s}}^{\prime }$ values that are considerably different between glossy and rough surface. This finding is again stressed by the red markers denoted as “no $R_S$” in Fig. 7(c), which displays box plots of the outlined glossy region (A′) and the rough region (B′). Besides the systematically larger values compared to the reference (dashed gray line), a difference in ${\mu }_{\mathrm{s}}^{\prime }$ of nearly $0.25{\mathrm{mm}}^{-1}$ between the two regions can be made out. By contrast, the evaluation with $R_S$, displayed in Fig. 7(a), only shows artifacts at the boundaries but the reduced scattering coefficients of glossy and rough surface areas are almost the same and agree with the reference values. These circumstances are again highlighted with boxplots displayed in Fig. 7(c). The green markers denoted at the top of the graph as “with $R_S$” indicate the averaged reduced scattering coefficient over the glossy region (A) and the rough area (B), both lie within the trust region of the known reference.

Fig. 7

(a) Reduced scattering map calculated with the modified theory taking the $R_S$ parameter into account. The labeled region (A) is glossy, whereas region (B) has a sand blasted rough surface. (b) Reduced scattering calculated with the state-of-the-art method without taking surface scattering into account, again the regions labeled (A′) and (B′) indicate the same different surfaces as in (a). (c) Box plot of the averaged reduced scattering coefficients of the regions (A) and (B), where the surface scattering model was utilized compared to (A′) and (B′) for which the surface scattering was neglected.

The corresponding maps and box plots are shown for ${\mu }_{\mathrm{a}}$ in Fig. 8. The contrast between glossy and rough surface for both evaluation methods “with $R_S$” and “no $R_S$” [compare Figs. 8(a) and 8(b)] is not as distinct as for ${\mu }_{\mathrm{s}}^{\prime }$ but still present. The median absorption coefficients ${\mu }_{\mathrm{a}}$ of the regions (A), (B), and (A′), (B′) charted in the boxplot [Fig. 8(c)] behave similar to ${\mu }_{\mathrm{s}}^{\prime }$. Without $R_S$, the absorption in general is systematically too high and the difference between glossy (A′) and rough (B′) surface is approximately $\approx 0.001{\mathrm{mm}}^{-1}$.

Fig. 8

(a) Absorption map calculated with the modified theory taking the $R_S$ parameter into account. The labeled region (A) is glossy, whereas region (B) has a sand blasted rough surface. (b) Absorption map calculated with the state-of-the-art method without taking surface scattering into account, again the regions labeled (A′) and (B′) indicate the same different surfaces as in panel (a). (c) Box plot of the averaged absorption coefficients of the regions A and B, where the surface scattering model was utilized compared to (A′) and (B′) for which surface scattering was neglected.

5.

Summary and Outlook

A modified solution of the RTE, which takes surface scattering into account, was presented. A parameter $R_S$ for quantifying surface scattering was introduced and the influence of this so-called surface scattering parameter on the determination of the optical properties was discussed. The modified theory was validated with different kinds of phantoms using SFDI experiments. These results, presented in Sec. 4, demonstrate significant improvements in the correctness and robustness of solving the inverse problem with the new model. Thus, a correction for systematic errors originated by diffuse surface reflections could be achieved without the need for crossed polarizers and consequently, a polarized light propagation theory.⁵ Since the surface roughness of the produced phantoms is comparable to that of skin, it is also demonstrated that surface roughness should be taken into account when applying SFDI to skin, e.g., for diagnosis.

In addition to the improved determination of the optical properties, the surface scattering parameter $R_S$ allows a quantification of the surface scattering. This could, for example, be used to differentiate between a lesion and surrounding healthy tissue, as the surface roughness in general is different.²⁶

Currently, the Lambert approximation does not allow for quantifying the surface roughness with respect to parameters usually applied in metrology. In future work, a more precise surface BRDF model should be applied in order not only to determine a correction parameter but, for example, the RMS surface roughness $R_q$, an important parameter in surface metrology. This would allow the calculation of ISO roughness maps. However, for such evaluations, more experimental data are needed. This is why the setup has to be modified such that angularly resolved measurements at least at two different angles become possible.