1. Introduction

To accurately model how photons propagate in tissue-like optical phantoms consisting of scatterers, absorbers, and bulk media [1], it is essential to gain knowledge on the optical properties, i.e. the refractive index, the scattering and absorption coefficients (μ _s and μ _a), and the scattering phase function (p(θ)). At distances of several d _mfp′ =1(μ _a+μ _s′) away from the source, where diffusion models are valid, only the reduced scattering coefficient μ _s′ = μ _s (1-g) is needed to describe light scattering [2]. For an accurate description close to the source using Monte Carlo (MC) simulations, the anisotropy factor g = g ₁ = 〈cos(θ)〉 (first moment) and higher order of moments (g_n,n≥2) of the phase function need to be taken into account [3–5].

Polystyrene microsphere suspensions are frequently used as the scattering component in optical phantoms. When using microspheres, Mie theory can be employed to accurately describe both the scattering coefficient and the phase function of a monodispersive suspension [6] or combinations of spheres with different, known sizes [4]. Michels et al. [7] used Mie theory, a sophisticated goniometry setup, spatially resolved reflectance, and collimated transmission to estimate the scattering properties of six different nutritional fat emulsions including their scattering size distribution of the scattering particles. However, in the latter study, the scatterer size distribution could not be freely fitted but had to be restricted to a linear decay in the logarithmic size scale.

The scattering properties of milk depend on the size distribution of fat globules and casein micelles, which varies greatly between homogenized and unhomogenized milk [8], and the refractive index differences between fat globules and casein micelles [9]. Hence, Mie theory cannot readily be employed to describe its scattering properties. Thueler et al. [4] proposed an adaptive light scattering model that includes a parameter γ = (1 - g ₂)/ (1 - g ₁) taking into account not only the anisotropy factor, but also the second moment of the phase function. Estimations of the scattering properties were achieved by matching the model with measurements from a fiber optic probe with multiple fiber separations (0.3 – 1.35 mm). Their model did not match the measurements for the closest two fibers (< 0.4 mm distance), when measuring on phantoms consisting of combinations of five polystyrene spheres, ranging from 0.05 μm to 1 μm in diameter. This indicates the need to include higher order of moments of the phase function to accurately model how photons propagate in microsphere suspensions.

The aim of this study was to present a method for determining the scattering properties of non-absorbing polydispersive optical phantoms in sufficient detail to allow accurate modeling of light intensities close to the source. The basis of the proposed method has some similarities with those presented by Thueler et al. [4], Palmer et al. [6], and Michels et al [7], such as: 1) The phase function estimation is improved by first measuring the scattering coefficient μ _s with a spectral collimated transmission (SCT) setup; 2) Measured and simulated spatially resolved diffuse reflectance (SRDR) with small source-detector separations are matched using non-linear fitting.

On the other hand, the proposed method differs from the methods presented in the above mentioned studies, in two important ways: 1) Mie theory is not used but the phase function is approximated with a two-parametric (α _Gk and g, see below) Gegenbauer kernel (Gk) phase function [10], enabling application to various polydispersive non-absorbing media with unknown size distributions, covering the γ range used by Thueler et al [4]; 2) Intensity calibration (simulations vs. measurements) is done either using a polystyrene sphere suspension (absolute calibration; no fitting) or in the fitting procedure (relative calibration). The latter does not, in contrast to Thueler et al. [4] and Palmer et al. [6], require measurements on reference microsphere suspensions. The robustness is further increased by measurements on multiple phantoms with different concentrations of scatterers. Together, this improves the accuracy in the two-parametric phase function determination in a way that enables accurate light modeling close to the source.

2. Material and methods

2.1 SCT setup

Collimated transmission is a widely used and accepted method for measuring the scattering coefficient (μ _s), absorption coefficient (μ _a) or the total attenuation coefficient (μ _t) of optical phantoms. The method is based on the Beer-Lambert law which describes collimated light attenuation as a function of optical pathlength through a medium. The essentials of the method can be found elsewhere [11]. A schematic view of the SCT setup is shown in Fig. 1 left panel. The two fiber collimators (Edmund Optics Inc., UK) were connected to a broadband tungsten halogen lamp (HL-2000, Ocean Optics Inc., USA) and a spectrometer (400–900 nm AvaSpec 2048-5-RM, Avantes BV, The Netherlands, grating: VB 600 lines/mm), respectively, using glass fibers with core/cladding sizes of 200/230 μm and a numerical aperture (NA) of 0.37. The two collimators, aligned and placed 250 mm apart facing each other, ensured a uniform and collimated illumination and detection of the light passing through the sample cuvette. The divergence of the collimators was 0.37 degrees. For each sample, 10 consecutive measurements where conducted at different sample thicknesses, achieved using a submerged glass cylinder connected to a high resolution translational micro stage. For each sample, the step size was adjusted to ensure that no measurement was conducted on optical thicknesses greater than 3 mfp (mean free path; mfp = 1/(μ _a+μ _s)).

2.2 SRDR setup

For the SRDR measurements, the radial decay in backscattered light intensity was measured using a fiber probe with several light collecting fibers positioned at different distances away from a light emitting fiber. The SRDR profile (magnitude and spatial shape) depends on the optical properties and the inverse problem can be solved using either the diffusion approximation [12] or Monte Carlo simulations [4].

A schematic view of the SRDR probe tip is shown in Fig. 1, right panel. The fiber-optic probe used in the SRDR setup consisted of 5 linearly aligned fibers, with one of the outermost fibers as the illuminating source fiber (S in Fig. 1 right panel) and the remaining 4 fibers as detector fibers (D _1–4 in Fig. 1 right panel). The center-to-center source-detector separations were 230 μm, 690 μ, 1150 μm, and 1610 μm, respectively. The spectrometer, light source and optical fiber specifications were the same as used in the SCT setup. The dependency between emission angle and light intensity of the emitting fiber was quantified in a goniometric setup. This angular dependency was found to correspond well with the specified NA, and was taken into account in the analysis of the simulation data by weighting each photon according to its emission angle. Similarly, the angle sensitivity of the backscattered light, collected by the fiber-optic probe when attached to the spectroscope, was characterized and taken into account in the analysis of the simulation data. This sensitivity, however, was found limited by the spectroscope to about ±6 degrees. When performing the measurements, the SRDR probe was handheld and vertically submerged (approx. 1 cm) in a 100 mL beaker. No movement artifacts due to inhomogeneous illumination were observed during the measurements.

 

Fig. 1. Left Panel: The SCT setup. Right Panel: The tip of the probe used for SRDR measurements, with source fiber (S) and detector fibers (D _1–4). Note: The SRDR probe is not drawn to scale, except for fiber separations. The distance from the edge of the probe to the edge of the closest fiber (S) was 2 mm.

Download Full Size | PDF

2.3 Spectral data preprocessing

Since the sensitivity of the spectrometer was insufficient below 500 nm and above 800 nm, only the 500–800 nm range was used in the analysis. For both the SCT and the SRDR measurements, dark spectra for every channel were recorded after each measurement and subtracted from the signal spectra before further processing. To compensate for any changes in lamp intensity and color during the SRDR measurements, all signal spectra were normalized with simultaneously recorded lamp spectra. Additionally, all measurements were normalized with a lamp-normalized white reference spectra recorded on a smooth surface covered with BaSO₄ with a reflectivity higher than 98% in the 500-800 nm range [13].

2.4 Preparation and scattering coefficient characterization of milk dilutions

Homogenized ultra-high temperature (UHT) treated milk (Arla Foods, Sweden), with a fat content of 1.5% and protein content of 3.4%, was used as a scattering basis for the optical phantoms. By diluting the milk with water, using a precision scale and assuming a density of 1.0 g/cm³ for the phantom components, four aqueous dilutions with different fractional milk contents were prepared. The dilutions were named OP ₁ through OP ₄ (OP = Optical Phantom), and the fractional milk contents were set to 1/8, 1/4, 1/2 and 1, respectively. The spectral scattering coefficient was determined by SCT measurements, assuming,μ _a = 0.0 mm^-1. The water absorption, with a maximum μ _a of below 0.003 mm^-1 in the wavelength region of interest [14], was found to have minimal or no impact on the detected intensity and was thus neglected.

2.5 Calibration of the SRDR system

To individually calibrate the detector fibers of the SRDR system, simulations and measurements of the diffuse reflectance spectra for an aqueous monodispersive polystyrene microsphere suspension (Bangs Laboratories, Inc., USA, radius 155 μm) were compared. The Monte Carlo (MC) setup consisted of a semi-infinite slab representing the suspension, with,μ _a = 0.0 mm^-1, μ _s as determined from a SCT measurement, and a Mie phase function calculated using a wavelength dependent relative refractive index based on water [15] and polystyrene [16]. The absolute calibration factors (A _abs) for each detecting fiber were calculated as the ratio between the mean intensity in the wavelength range 500–800 nm of the measured and the simulated spectra.

The relative calibration factors (A _rel) between the detecting fibers were determined from a measurement where all fibers were illuminated uniformly and equally. For this measurement, the probe was vertically submerged in a 100 mL beaker containing a high-scattering polystyrene sphere suspension, which was illuminated from below by a white halogen lamp (Flexilux 150 HL, Schölly Fiberoptik GmbH, Denzlingen, Germany).

2.6 Two-parametric phase function determination

The idea for determining the two-parametric phase function was to solve the inverse problem by comparing measured SRDR spectra with pre-simulated MC spectra. The MC setup used the two parametric Gk phase function (parameters α _Gk and α _Gk):

The anisotropy factor g depends on both α _Gk and g _Gk; g _Gk only equals g when α _Gk = 0.5 (i.e. the Henyey-Greenstein phase function) [10]. Simulations were performed for 7 values of α _Gk (0.025, 0.05, 0.1, 0.2, 0.3, 0.4, and 0.5) and 11 values of the anisotropy factor g (0.7, 0.72, …, 0.9). The g _Gk parameter needed for simulation input was numerically calculated using the relationship between g _Gk, α _Gk and g described in [10]. Each of these 7 × 11 = 77 simulations were performed by injecting photons in a single point. The direction of each photon was randomly chosen based on a uniformly illuminated solid angle (limited by NA = 0.37) assumption. Only a single scattering coefficient, denoted μ _s,sim =3.72 mm^-1 (average μ _s for OP ₂ in the wavelength range 500–800 nm), was simulated. Additional μ _s levels were added in the post-processing by rescaling the original simulation data [17]. In short, the rescaling is performed by multiplying the optical pathlength and the photon detection point from the original simulation by the factor μ _ssim/μ _s . For each of the 77 combinations of α _Gk and g, recalculations were done for the SCT measured μ _s levels at wavelengths λ = 500, 510, …, 800 nm for each of the four phantoms (OP _1–4). The effect of the source distribution (radius 100 μm) was added after rescaling by randomly offsetting each detection position. The emission and detection angle dependency of the SRDR system was taken into account by adjusting the weight of each photon accordingly, as described earlier. Finally, the detected intensity at the 4 source-detector separations (D _1–4) was calculated, based on at least 1.0×10⁵ detected photons.

Altogether, this resulted in a MC generated multi-dimensional intensity array with the dimensions (lengths) α _Gk(7), g(11), λ (31), OP _1–4(4), and D _1–4(4). Finally, the intensity array was smoothed in the α _Gk and g dimensions using cubic polynomials to eliminate local minima caused by small stochastic variations in the simulation data. No trends in the residuals between smoothed and original data could be seen, and the maximum deviation was below 2%. The final intensity array was denoted I _MC. Similarly, the array containing the measured intensities was denoted I _meas, with the dimensions (lengths) λ(31), OP _1–4 (4), and D _1–4 (4).

The error functions between MC generated and measured intensity (Eq. 2: absolute calibration method and Eq. 3: relative calibration method) were minimized in a least-square sense using the Levenberg-Marquardt algorithm with α_Gk and g as fitting parameters. The initial values of these parameters were chosen as the mid-points in the two vectors α _Gk and g (0.2625 and 0.8), respectively. Evaluating other initial values yielded identical results. Both α _Gk and g are assumed to depend only on λ, and not OP and D .

The minimizers were denoted α _{Gk, abs} ^* (α _{Gk, rel} ^*) and g ^* _abs(g ^{* rel}), giving the optimal MC generated intensity, I ^* _MC,abs (I ^* _MC,rel), for the absolute (relative) calibrated case. Note that the minimizers are wavelength dependent and that the optimization based on the relative error function needs to be performed for all wavelengths simultaneously.

For the phase function approximation, OP ₄ was omitted from the analysis, since a high concentration of scatterers likely affects the phase function (see Discussion). In addition, for some optical phantoms and wavelengths, the signal level for D ₄ was below two standard deviations of the detector dark noise and was therefore excluded from the analysis.

The robustness of the method was evaluated by excluding data from the optimization procedure. This was achieved by limiting the measured data utilized for the phase function approximation to either a single fiber separation ([OP _1–3, D_j], j = 1, 2, 3) or a single optical phantom ([OP_i, D _1–3], i = 1, 2, 3), resulting in six different α _Gk,abs ^* and g ^* _abs combinations, and six different α _Gk,rel ^* and g ^* _abscombinations. For each of these six combinations, I ^* _MC,abs and I ^* _MC,rel were evaluated for all nine OP_i and D_j combinations and compared to the corresponding I _meas. The deviation between simulated and measured spectra was quantified as the root-mean-square (rms) errors

and

For comparison, the optimal anisotropy factor g ^* _HG,abs for a Henyey-Greenstein (HG) phase function approximation (α _Gk = 0.5) was calculated using [OP _1–3, D _1–3].

2.7 Optical phantoms including absorption

To investigate how the two-parametric phase function describe light propagation in optical phantoms containing both absorption and scattering components, SRDR measurements were performed on three scattering and absorbing optical phantoms (OP _5–7), and were compared to simulated spectra. As absorption component, water soluble, oil-free blue ink (Coloris, Germany) was used. The ink was diluted with water to 1/500 of the original concentration, and the spectral absorption coefficient was determined by SCT measurements. The fractional content of the diluted ink in OP _5–7 was set to 1/20, 1/7 and 1/2, respectively. The fractional milk concentration was set to 1/4 (equal to the scattering component concentration of OP ₂) in all three optical phantoms.

The simulated spectra for OP _5–7 were calculated based on both (α _Gk,abs ^*, g ^* _abs) and (α _Gk,rel ^*, g ^* _rel), determined from [OP _1–3, D _1–3]. By logging the optical pathlength of every detected photon in the original simulation, inclusion of the absorption effect was done after the μ _s,sim/μ _s rescaling by reweighting every detected photon by a factor of exp(-μ _a·PL), where PL is the rescaled optical pathlength [17]. The μ _a values used in the simulated spectra were obtained by multiplying the fractional content described above by the μ _a spectra from the SCT measurements of the blue ink dilution.

3. Results

3.1 Calibration of the SRDR system

The fit between the measured and simulated spectra for the polystyrene microsphere suspension is shown in Fig. 2 left panel. The measured spectra are multiplied with the absolute calibration factors A _abs. The right panel shows a comparison between A _abs and A _rel. In this figure, both factors were normalized with its mean over D. The quotients (A _rel/〈A _rel〉_D)/(A _abs/(〈A _abs〉_D) were -2.2%, -0.6%, 1.3%, and 1.0% for D ₁, D ₂, D ₃ and D ₄, respectively.

 

Fig. 2. Left Panel: Calibration measurements using a polystyrene microsphere suspension. Measured (black) and simulated (gray) spectra for the source-detector separations D _1–4. Measured spectra are normalized by A _abs (see Eq. 2). Note: The intensities are unitless, since the simulated spectra are normalized with the incident light intensity. Right Panel: A _abs for D _1–4 from the microsphere suspension measurement (*) and A _rel from the uniform illumination measurement (Δ). Both series are normalized by their mean value.

Download Full Size | PDF

3.2 Scattering coefficient

The maximal optical path lengths during the measurements were in the interval 0.15 mm (OP ₄) and 0.55 mm (OP ₁), giving an mfp < 3 (for OP ₃ and λ =500 nm). The small acceptance angle of the fiber collimator (0.37 degrees), makes detection of scattered light negligible. Calculations according to Wang et al. 11 or OP ₃ and the Gk phase function with α _Gk,abs ^* and g ^* _abs (Fig. 4), gives a μ _s -error less than 0.3%. In Fig. 3 left panel, μ _s, as measured by SCT, is shown as a function of wavelength for OP _1–4. The scattering coefficient of the optical phantoms displayed a linear response to concentration, except for OP ₄, which displayed lower values, as seen in the quotients μ_s,OPi μ_s,OP1 OP 1, i = 2, 3, 4, Fig. 3 right panel. This non-linear response between μ _s vs. concentration has been observed by others [18–20], and is further commented in the discussion.

 

Fig. 3. Left Panel: μ _s for OP _1–4 as a function of wavelength. Right Panel: μ_s,OPi/μ_s,OP1, i = 2,3,4 . Dotted lines: Expected quotients assuming linear response between μ _s and concentration. Solid lines: experimental quotients.

Download Full Size | PDF

3.3 Two-parametric phase function determination

In Fig. 4 left panel, α ^*_Gk,abs and α^* _Gk,rel, utilizing [OP _1–3, D _1–3], are shown, and g ^* _abs and g ^* _rel are displayed in Fig. 4 right panel.

 

Fig. 4. The wavelength-resolved parameters for the two-parametric phase function. Left Panel: α _Gk,abs ^*, () and α _Gk,rel ^* (). Right Panel: g ^* _abs (), g ^* _rel() and g _HG,abs (O).

Download Full Size | PDF

For both the absolute calibrated and the relative calibrated case, the anisotropy factor displayed decreasing values with increasing wavelength. The anisotropy factor with the HG phase function was lower than those with the Gk phase function. The difference between g ^* _rel and g ^* _abs was in the interval -1.2% to -0.8%.

The range and the average of the rms-errors, as given by Eq. 4 and 5, for each of the six Gk phase functions estimated from the subsets [OP _1–3, D _j], j = 1, 2, 3, and [OP_i, D _1–3], i = 1, 2, 3, are presented in table 1.

Table 1:. Deviations (E ^rms _abs and E ^rms _rel) between simulated and measured spectra evaluated over OP _1–3, D _1–3, using different data-sets ^a) in the phase function estimation. Deviations are given as mean and min-max range.

View Table

When utilizing all available data [OP _1–3, D _1–3] in the phase function approximation, the worst case rms-errors in Eq. 4 and 5, were 6.5% and 2.5% for the absolute and relative calibrated cases, respectively. For the absolute calibrated case, the corresponding value when utilizing the HG phase function and [OP _1–3, D _1–3], were 8.5%, with a min - max range of 3.3% to 16.8%, displaying substantially larger errors compared to the cases utilizing the Gk phase function.

For phase function approximations utilizing only subsets of OP and D, Eq. 4 displayed worst case rms-errors above 10% except when calibrating using [OP ₂, D _1–3]. In general, the relative calibration case displayed lower rms-errors, even when the determination was based on a sub-set of data while evaluation was performed on the complete data-set. Fig. 5 shows measured intensity, I _meas, together with I _MC,abs ^* and I _MC,rel ^* based on the phase function approximation utilizing [OP _1–3, D _1–3]. There was no relationship between the rms-error and either of the level of scattering (OP_i, i = 1, 2, 3) or the source-detector distances.

 

Fig. 5: Measured (black) and simulated (solid gray : I _MC,abs ^*, dashed black : I _MC,rel ^*) intensities, with phase function approximation based on [OP _1–3, D _1–3]. Left Panel: D ₁, Right Panel: D ₃.

Download Full Size | PDF

3.4 Optical phantoms including absorption

In Fig. 6 left panel, the absorption coefficient of the diluted (1/500) blue ink as a function of wavelength is shown. Fig. 6 right panel displays the comparison between measured and simulated spectra in D ₃ for OP ₅ and OP ₇. D ₃ was chosen as an illustrative example in favor of D ₁ and D ₂, since the effect of absorption is most noticeable in D ₃.

 

Fig. 6: Left Panel: μ _a for diluted blue ink (1/500) as function of wavelength. Right Panel: Measured (black) and simulated (solid gray : I _MC,abs ^*, dashed black : I _MC,rel ^*) intensities for D ₃.

Download Full Size | PDF

For OP _5–7 the rms-errors were below 13.4% and 7.3% for the absolute and relative calibrated cases, respectively. There was no relationship between the rms-error and the source-detector distance, while the error increased with the level of absorption (in OP_i, i = 5, 6, 7). The only rms-errors above 10% were for OP ₇ and the absolute calibrated cases.

4. Discussion and conclusions

We have presented a method for a two-parametric phase function determination at small source-detector separations that can be applied to polydispersive scattering liquids. The method involves spectral scattering coefficient determination, SRDR measurements, and MC modeling and simulation of light transport. The phase function parameters are determined by nonlinear fitting minimizing the difference between measured and simulated SRDR data. To speed-up the inverse Monte Carlo algorithm, pre-simulated multi-dimensional data-sets, based on a two-parametric phase function, and scaling for μ _s, is utilized. The inverse algorithm employs spectral scattering coefficients that are separately determined using collimated transmission measurements.

Evaluations were performed on four different concentrations of UHT-treated milk with 1.5% fat content and 3.4% protein content. Two different intensity calibration methods were presented; one using a monodispersive polystyrene sphere suspension for an absolute calibration, and another using a relative calibration scheme where the difference in sensitivity between fibers is equalized.

4.1 Spectral scattering coefficient

The measured spectral scattering coefficients for OP _1–3 were found to be multiples of each other which is expected since the fractional milk content in OP _1–4 was set to 1/8, 1/4, 1/2, and 1, respectively. Assuming a linear relationship between concentration and μ _s, the measured μ _s for OP ₄ was smaller than expected. However, several studies [19, 20] have reported that for sufficiently high volume fraction of scatterers in liquid phantoms, the scattering processes are not independent of each other. Zaccanti et al. [20] show that deviations in μ _s are larger than for μ _s′, indicating changes in both μ _s and the anisotropy factor. In their work, this nonlinear response was observed for volume fractions larger than 0.01. In our work, the volume fraction of scattering particles (fat globules and casein micelles) in OP ₄ (undiluted milk) was calculated to be 0.042, assuming 80% of the protein to be casein micelles [21]. The milk concentration for OP ₄ clearly was not in this independent-scattering concentration region, which was the reason for exclusion in the phase function approximation. For OP _1–3 the measured μ _s displayed a mean (over wavelength) deviation of less than 1 % from the expected μ _s when assuming independent scattering events.

4.2 Calibration of the SRDR system

The proposed method involves intensity calibration on a single known optical phantom (absolute calibration method) or the calibration is done in the fitting procedure (relative calibration method), in contrast to several other calibrated SRDR methods [22–26], which relies on an extensive number of phantoms. Most of those methods are also based on adaptations of analytical light transport models, which have shown to fail when dealing with small source-detector separations or large ranges in optical properties.

In the absolute calibration method, an alteration in the calibration factors will affect the absolute level of both α _Gk,abs ^* and g ^* _abs. However, repeated calibration measurements displayed intensity variations not greater than 5%. Tests have shown that a corresponding alteration in the intensity of the measured spectra gave a change in calculated anisotropy within the interval 0.005 to 0.01 (varying somewhat with wavelength), indicating a robust calibration process.

In the relative calibration method, the amplification factors A _rel only compensated for inter-channel differences in fiber transmittance, fiber-spectroscope coupling, and detector efficiencies. The scaling factor between the measured and MC simulated SRDR intensities is compensated for by the normalization described in Eq. 3. As a measure of how well the relative calibration method performs compared to the absolute calibrated method, we calculated the scaling factor k between I _meas and I _MC,rel ^*. When comparing 〈k A _rel〉_OP,D to 〈A _abs〉_OP,D, the difference was -5.8%. This difference was only slightly larger than the intensity variations for the repeated calibration measurements in the absolute calibration method, suggesting that the relative calibration performed equally well as the absolute calibration method, given that [OP _1–3, D _1–3] is used.

4.3 Two-parametric phase function determination

The Gk phase function was determined utilizing seven different combinations of OP and D, data for both the absolute and relative calibrated method. In addition, the anisotropy factor was calculated based on a HG phase function approximation. Other authors [24, 27] have used Mie theory when describing light scattering in liquid optical phantoms. However, if the Mie phase function would be employed in the MC simulations used in our phase function approximation algorithm, every particle size would correspond to one dimension in I _MC. In liquids containing many different scattering particle sizes, such as milk, the amount of data in this array would be huge, and the overall MC simulation time becomes unrealistic. This could be solved by using a semi-analytic Monte Carlo approach where the simulation data is rescaled in each optimization iteration [28]. However, it is likely that multiple solutions to the inverse problems are rapidly introduced as the number of particle sizes increases. The difficulty of finding a scatterer size distribution has been observed by others [7].

When utilizing [OP _1–3, D _1–3] for the Gk phase function determination, the results for α _Gk and g were very similar in both spectral behavior and absolute levels when comparing the absolute and relative calibrated methods. The I _MC ^* vs. I _meas rms-errors were lower using the relative calibrated method, except when utilizing [OP ₃, D _1–3] for the phase function approximation. For phase function approximations utilizing only subsets of OP and D, the I _MC ^* vs. I _meas rms-errors were larger compared to the case when utilizing [OP _1–3, D _1–3]. In general, utilization of only D ₁ in the phase function approximation resulted in I _MC ^* < I _meas, whereas I _MC ^* > I _meas when only D ₃ was utilized for the phase function approximation. In addition, the D ₃ intensity variation in I _MC as a function of α _Gk was very small. For example, inspection of I _MC(α _Gk, g = 0.80, λ = 500nm, OP ₁, D) showed an intensity variation of 31 % for D ₁, 4.6 % for D ₂ and 3.3 % for D ₃. It is, therefore, reasonable to assume that for the scattering components and the concentrations used in this study, determination of the two-parametric phase function has to be performed at small source-detector separations (0.23 mm) and not larger than 1.2 mm.

The comparison between simulated (both absolute and relative calibrated) spectra and measured spectra for the optical phantoms including absorption (OP _5–7), gave only slightly larger average rms-errors than in the OP _1-3 cases. These results are an additional indicator that the two-parametric phase function as determined by the proposed method is able to predict the light propagation also in liquid optical phantoms including absorption. The absorption level in OP ₇ is close to 1 mm^-1, which is, in terms of biological tissue, regarded as very high.

Anisotropy values reported by other authors range from 0.70 to 0.86 (633 nm to 700 nm), using techniques such as oblique angle illumination [29], coherent backscattering measurements [30], or integrating sphere measurements [31]. This is in agreement with our findings, where the anisotropy factor ranged from 0.79 to 0.86, utilizing the Gk phase function, and 0.75 to 0.82, utilizing the HG phase function. However, it should be noted that the values varies with temperature, fat content, different pasteurization and homogenization processes, dilutions etc. Therefore, it is recommended that when using milk or fat emulsions [7] as a scattering component in optical phantoms, the proposed procedure should be performed for the specific batch of scatterers at hand.

In conclusion, we have presented a method for estimating a phase function that accurately describes how light propagates in polydispersive scattering liquids. A multi-dimensional Monte Carlo model based on a two parametric Gegenbauer-kernel phase function approximation and μ _s values derived from SCT, was used to simulate a multi-dimensional array of SRDR spectra. The two phase function parameters were estimated from minimizing the rms-error between simulated and measured SRDR spectra using a non-linear search algorithm. Intensity calibration of the SRDR measurements was done either by a separate measurement on a monodispersive polystyrene sphere aqueous suspension (compensating for measured-simulated amplification differences) or by uniformly illuminating the fiber probe (only compensating for between-fiber amplification differences). An additional overall amplification factor, calculated as the average ratio between fitted and measured spectra, was added to the relative calibration case when comparing intensities. Both the absolute and the relative calibrations yielded similar phase function estimations when measuring on several non-absorbing phantoms containing different concentration levels of the same scattering substance. The deviation between fitted and measured spectra was low, even when adding ink to the samples, which strongly supports that the estimated phase functions are accurate.