Spatial frequency domain imaging using an analytical model for separation of surface and volume scattering

Abstract. A method to correct for surface scattering in spatial frequency domain imaging (SFDI) is presented. The use of a modified analytical solution of the radiative transfer equation allows calculation of the reflectance and the phase of a rough semi-infinite geometry so that both spatial frequency domain reflectance and phase can be applied for precise retrieval of the bulk optical properties and the surface scattering. For validation of the method, phantoms with different surface roughness were produced. Contrarily, with the modified theory, it was possible to dramatically reduce systematic errors due to surface scattering. The evaluation of these measurements with the state-of-the-art theory and measuring modality, i.e., using crossed linear polarizers, reveals large errors in the determined optical properties, depending on the surface roughness, of up to ≈100%. These results were confirmed with SFDI measurements on a phantom that has a structured rough surface.


Introduction
Quantitative spatial frequency domain imaging (SFDI) for separation of the absorption (μ a ) and the reduced scattering coefficient (μ 0 s ) 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][2][3] To eliminate specular and surface reflections in SFDI, crossed linear polarizing filters are generally used. 4 However, as shown by Wiest et al., 5 the use of crossed polarizers causes significant errors in μ 0 s and μ 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 SFD ðFÞ. This offset mostly affects high spatial frequencies F, as for those R SFD ðFÞ is steadily decreasing 6,7 and therefore, the ratio R s R SFD ðFÞ is increasing. Correcting for this offset without any more information than R SFD ðFÞ is hardly possible, because the subdiffusive scattering properties (e.g., phase function) and the surface scattering affect the spatialfrequency resolved reflectance in a similar way. Bassi et al. 8 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. 9 First, results associated to this work were already presented on the European Conference on Biomedical Optics in 2017. 10 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 SFD ðFÞ but also the spatial shift or phase ϕ 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.

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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 2 7 0 R med ðq; n med ; n ext ; μ a ; μ 0 s ; ρ V ; θ 0 0 ; θ NA Þ; (1) originating from the volume scattered light. Note that the solution in Ref. 12 has its source inside the medium, this is why in case of an external source, it additionally has to be scaled by a transmission coefficient T V . The reflectance depends on the angular spatial frequency q ¼ 2πF, the refractive indices of both the internal (n med ) and external (n ext ) medium, the absorption coefficient μ a , the scattering μ s , or alternatively, the reduced scattering coefficient μ 0 s , the phase function ρ V of the scatterers and the angle: which is the angle inside the medium calculated according to Snell's law from the angle θ 0 under which the light is irradiated relative to the normal of the sample surface. Moreover, the limited detection angle with respect to the numerical aperture of the objective lens is taken into account with NA ¼ sinðθ NA Þ.
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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 4 4 5 I I ðx; y; zÞ ¼ with different spatial frequencies q i . The additional phase ϕ S ðq; zÞ is due to the topography of the sample, which in detail is explained later in this section. The reflected intensity I V;q i ðx; y; zÞ due to volume scattering then yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; 3 6 0 where the transmission coefficient E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 ; 6 3 ; 2 6 6 T V ≡ 1 −R S is the fraction of light that is not reflected at the surface but transmitted into the medium. The total surface reflectionR S is defined as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; 2 0 6R E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 1 6 1 which is the ratio of the integral over all directions of the surface reflected radiance L S and the integral over the half-space Ω of the incident radiance L I . The transmitted light propagates through the turbid medium according to the RTE, whereby the modulation amplitude [thin red line in Fig. 1(a)] of the light being reflected from the sample can be expressed in terms of the analytical solution Eq. (1) of the RTE as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 5 0 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 3 2 6 ; 4 5 5 In this equation, the radiance L V is, as illustrated in Fig. 1(b), originated by light that was scattered in the turbid volume. However, besides the often investigated change in the amplitude R V , scattering also leads to a significant spatial shift x ϕ ðqÞ of the detected sinusoidal intensity compared to the irradiated pattern. 8 This phase shift shall be discussed shortly in order to understand the principle of separating volume and surface scattered light, which will be presented afterward. The reason for this shift is that the reflected light in average propagates through the medium by the effective free path before the first interaction with the medium takes place. The oblique irradiation then results in a projection of this mean free path to the x dimension, which corresponds to the measured spatial shift x ϕ ðqÞ. The phase E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 2 6 2 ϕ med ðqÞ ¼ arg½R med ðq; n med ; n ext ; μ a ; μ 0 s ; ρ; θ 0 ; θ NA Þ; (9) due to the volume scattering therefore depends on both the incident angle θ 0 and the optical properties. However, experimentally, the phase shift always refers to a reference measurement, whereas the surface of the reference sample defines the position z ¼ 0. Any type of misalignment z ≠ 0 of the sample surface then will cause an additional phase shift of E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 1 6 5 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 (a) (b) Fig. 1 (a) Schematic illustration of the radiance inside a turbid medium, which is obliquely illuminated at θ 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.
Journal of Biomedical Optics 071604-2 July 2019 • Vol. 24 (7) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 7 5 2 x ϕ ¼ ϕ V ðqÞ q ; (11) with the definition: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 7 0 6 ϕ V ðq i Þ ≡ ½ϕ S ðq ¼ q i Þ þ ϕ med ðq ¼ q i Þ: (12) 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 Scratch-Dig < 80, 13 the total reflectionR S is given by Fresnel's equation according to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 5 9 2R S;Fresnel ðθ 0 ; n med ; n ext Þ ¼ with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 4 9 3 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 directioñ sðθ S ; ϕ S Þ with spherical coordinates θ S and ϕ S , the so-called bidirectional reflectance distribution function (BRDF): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 6 3 ; 3 3 3 ρ BRDF;S ðθ I ; ϕ I ; θ S ; ϕ S Þ ¼ dL S ðθ S ; ϕ S Þ L I ðθ I ; ϕ I Þ cosðθ i Þdθ I dϕ I (15) has to be introduced. The index S indicates that the light is only scattered by the surface and did not penetrate into the volume. This function describes the ratio between the incident radiance L I ðθ I ; ϕ I Þ and the differential radiance dL S ðθ S ; ϕ S Þ of surface scattered light being detected, where the incident directioñ sðθ I ; ϕ I Þ is given by the polar angle θ I and the azimuthal angle ϕ I . 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 14 showed such Lambert surfaces do not exist in practice, thus, more precise models have to be applied. For instance, Sun 15 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 9 6I S;q i ðx; y; zÞ ¼ R S I 0 2 ½sinðq i x þ ϕ S ðq i ; zÞÞ þ 1; (16) due to surface scattered light only depends on the surface scattering correction parameter R S . This parameter is according to Fig. 1

(b) given by the integral:
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 3 2 6 ; 7 1 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 6 6 0 and therefore, the exact knowledge of the BRDF is not necessary. The intensity defined in Eq. (16) therefore can be seen as this fraction of surface scattered light, which falls into the aperture of the detection system. In the case of a rough surface with a certain BRDF, the total reflection parameterR S is no longer given by Fresnel's equation for plane surfaces but instead reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 5 5 6R Within this work, the Lambert approximation ρ BRDF;S ¼ constant was used, for which Eq. (17) together with Eq. (19) yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 3 2 6 ; 4 8 0R 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∕μ 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., 16 when considering a stationary isotropic Gaussian surface, the KA can be applied, if

2.76λ
; where λ is the electromagnetic wavelength and k ¼ 2π∕λ the wave number.
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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 3 2 6 ; 9 2 I q i ðx; y; zÞ ¼ I V;q i þ I S;q i ; Journal of Biomedical Optics 071604-3 July 2019 • Vol. 24 (7) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 6 3 ; 7 4 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 6 3 ; 7 0 1 þR S I 0 2 ½sinðq i x þ ϕ S ðq i ; zÞÞ þ 1; (23) of these two detected sine shaped intensities I V;q i and I S;q i again is a sine (see derivation in Appendix A): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 6 3 ; 6 4 3 with the amplitude E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 6 3 ; and phase E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 6 3 ; 5 2 1 The amplitude measured for an angular spatial frequency q i , thus, yields I 0 2 R SFD ðq ¼ q i Þ, which in future is denoted as AC reflectance. The part I 0 2 ðR S þ R V ðq ¼ 0 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 ≠ 0 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 ϕ ¼ ϕ SFD ðq ¼ q i ; zÞ∕q i , where ϕ 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 SFD and phase ϕ SFD versus spatial frequency F ¼ q∕2π. The optical properties for this calculation were μ a ¼ 0.01 mm −1 , μ s ¼ 4.00 mm −1 , and n med ¼ 1.4, whereby a Henyey-Greenstein phase function with g ¼ 0.75 was utilized resulting in μ 0 s ¼ 1.00 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 SFD and the changes of the phase can clearly be seen.
In order to be able to neglect the influence of the phase function ρ V ðθÞ 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 ∈ ½0; M þ 1 expansion coefficients: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 3 2 6 ; 5 5 4 instead of the exact phase function. Though the accuracy of the solution depends on the expansion order, it was found that for an order of M ≥ 9, these deviations do not significantly influence the results. Thus, within this work, all solutions of the RTE were calculated with an order of M ¼ 9.

Materials and Methods
All SFDI measurements within this work were carried out with the spatial frequency domain setup described by Bodenschatz et al. 17   By projecting three sine patterns of J and different spatial frequencies F i , each phase shifted by 0, 2π∕3, and 4π∕3, images with intensity values I n;ðu;vÞ ðF i Þ for n ∈ ½1;2; 3 at pixel ðu; vÞ are acquired. Demodulating each pixel 18,19 yields the DC intensity as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 6 3 ; 6 9 7 I DC;ðu;vÞ ðF i Þ ¼ 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 mm −1 , Afterward, the model introduced in Sec. 2 was fitted to the data using a nonlinear least squares DOGLEG 20 algorithm, which was taken from the library CERES SOLVER. 21 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. 22 Nine different phantoms were produced, from which three are nonturbid with μ 0 s ≈ 0 mm −1 [ Fig. 3(a), pa0], the next three have reduced scattering coefficients of μ 0 s ≈ 1 mm −1 [ Fig. 3(a), pa1] and another three have μ 0 s ≈ 4 mm −1 [ Fig. 3(a), pa4]. Within one group of phantoms with same reduced (a) (b) 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 μ 0 s ≈ 0 mm to μ 0 s ≈ 4 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.
Journal of Biomedical Optics 071604-5 July 2019 • Vol. 24 (7) 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 mm length by scanning a diamond stylus (tip radius 5 μ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: According to the definitions given by the International Organization of Standardization (ISO), the so-called root mean square deviation R q and the maximum height of profile R Z were calculated from the surface profiles. These parameters are given in Table 1 for the three different surface roughnesses s1, s2, and s3. The values R q and R Z of the produced phantoms are within the range of values found in literature for skin. For example, Tchialeva et al. 23 gave a summary of the skin surface roughness obtained with different optical methods, which yield a R q value between 0.1 and 10 μm. 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 ½μ a ð600 nmÞ ≈ 10 −3 mm −1 , for which the absorption spectrum in detail is given by Krauter et al. 22 Additionally, one phantom with homogeneous optical properties of μ 0 s ≈ 1.05 mm −1 and μ a ≈ 0.44 mm −1 for the wavelength λ ¼ 700 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 AE 0.02Þ μm and R Z ¼ ð1.4 AE 0.3Þ μ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 AE 0.3Þ μm and R Z ¼ ð11.6 AE 0.5Þ μ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.

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.

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 SFD and ϕ SFD . The top plots in Table 1 The table gives the root mean square parameter R q and maximum height of profile R Z 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 l c , 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   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 ðmeasurement − fitÞ∕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 . Fitting the modified theory yields besides the optical properties μ 0 s and μ 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 SFD could be demodulated according to Eq. (29) and as R V ¼ 0 for these phantoms relation, Eq. (25) yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 ; 6 3 ; 4 9 9 R S ¼ R SFD : (33) These R S values were taken as reference and are denoted R S;ref . The trust region of AE3% 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;ref and the green squares indicate the results for the phantom group pa4, respectively. The dashed gray line indicates the target values with AE3% 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 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 μ 0 s ≈ 1 mm −1 and pa4 with μ 0 s ≈ 4 mm −1 . The plots in Fig. 5  display μ 0 s 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 AE5% for different wavelengths. The trust region of μ 0 s was estimated by the production accuracy of the epoxy resin phantoms according to Krauter et al. 22 . 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 AE0.25 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 μ 0 s 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.
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 μ 0 s 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 μ 0 s 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.
(b) (a) Fig. 5 The plots in (a) and (b) show the fitted reduced scattering coefficients μ 0 s 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 AE5% relative deviation.

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 ≈ 0.8% and for the rough surface at about R S ≈ 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 ≈ 5.5% for such surfaces, which also agrees to the R S values obtained for the rough sandblasted structure. 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 μ 0 s 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  Journal of Biomedical Optics 071604-8 July 2019 • Vol. 24 (7) region (B′). Besides the systematically larger values compared to the reference (dashed gray line), a difference in μ 0 s of nearly 0.25 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.
The corresponding maps and box plots are shown for μ 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 μ 0 s but still present. The median absorption coefficients μ a of the regions (A), (B), and (A′), (B′) charted in the boxplot [ Fig. 8(c)] behave similar to μ 0 s . Without R S , the absorption in general is systematically too high and the difference between glossy (A′) and rough (B′) surface is approximately ≈0.001 mm −1 .

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. 5 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. 26 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.

Appendix A: Sum of Two Sine Functions
It is to be shown that the sum of two sine functions F 1 ¼ R 1 sinðθ 1 Þ and F 2 ¼ R 2 sinðθ 2 Þ yields a sine function 27 of kind F 1;2 ¼ F 1 þ F 2 ¼ R 1;2 sinðθ 1;2 Þ with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 4 ; 3 2 6 ; 3 4 6 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 5 ; 3 2 6 ; 3 0 1 θ 1;2 ¼ arctan2ðR 1 cosðθ 1 Þ þ R 2 cosðθ 2 Þ; R 1 sinðθ 1 Þ þ R 2 sinðθ 2 ÞÞ: Consider R 1 , R 2 being the absolute values and θ 1 , θ 2 the arguments of the two complex numbers C 1 ¼ R 1 expðiθ 1 Þ and C 2 ¼ R 1 expðiθ 2 Þ, as shown in Fig. 9. The sum C 1;2 ¼ C 1 þ C 2 of these two complex numbers is just a sum of two vectors in the complex plane, for which the amount of C 1;2 is given by the law of cosine as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 6 ; 3 2 6 ; 1 9 0 The angle between the real axis and C 1;2 , which corresponds to the phase θ 1;2 , is coupled via the tangent with the imaginary: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 7 ; 3 2 6 ; 1 2 9 IðC 1;2 Þ ¼ IðC 1 Þ þ IðC 2 Þ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 8 ; 3 2 6 ; 9 9 ¼R 1 sinðθ 1 Þ þ R 2 sinðθ 2 Þ; and real part By applying the inverse of the tangent, the last expression yields as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 ; 6 3 ; 3 9 1 θ 1;2 ¼ arctan2ðR 1 cosðθ 1 Þ þ R 2 cosðθ 2 Þ; R 1 sinðθ 1 Þ þ R 2 sinðθ 2 ÞÞ: As the imaginary parts of C 1 , C 2 , and C 1;2 exactly correspond to the sine functions Eqs. (34) and (35) are proven.

Disclosures
The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.