Determination of three optical properties from subdiffusive spatially resolved reflectance calculations

Philipp Krauter; Dominik Reitzle; Simeon Geiger; Alwin Kienle

doi:10.1117/1.JBO.22.7.075003

21 July 2017 Determination of three optical properties from subdiffusive spatially resolved reflectance calculations

Philipp Krauter, Dominik Reitzle, Simeon Geiger, Alwin Kienle

Author Affiliations +

Funded by: Evangelisches Studienwerk Villigst

Journal of Biomedical Optics, Vol. 22, Issue 7, 075003 (July 2017). https://doi.org/10.1117/1.JBO.22.7.075003

Abstract

We report a theoretical study on the determination of three optical properties from spatially resolved reflectance calculations. In particular, the reduced scattering coefficient μ_s^′, the absorption coefficient μ_a, and the recently defined phase function parameter σ are identified. The solution of the inverse problem is based on the principal component analysis of a large set of reflectance profiles that were calculated using an analytical solution of the radiative transfer equation. Different phase function types were studied to test the method in the range of 0.63 mm⁻¹≤μ_s^′≤4.2 mm⁻¹ and 0.002 mm⁻¹≤μ_a≤0.1 mm⁻¹. For curves impaired with noise, we were able to reconstruct μ_s^′ and μ_a with relative median errors of 2.5% and 12%, respectively, and σ with an absolute median error of 0.04.

1. Introduction

Optical light is a powerful diagnostic tool for medicine and material analysis. Spatially resolved reflectance measurements are a common method to investigate light propagation in scattering media. The analysis can often be summarized statistically in terms of optical properties such as the reduced scattering coefficient $μ_{s}^{'}$ , the absorption coefficient $μ_{a}$ , the refractive index $n$ , and the scattering phase function $p (θ)$ . Although it has long been known that the phase function has significant influence on light propagation,¹^–³ it has often been neglected, in particular, for reflectance experiments. In most studies, either the underlying phase function was known and taken into account appropriately or an arbitrary phase function was assumed. Only in a few cases was a determination of a phase function parameter performed.²^,⁴^,⁵ In principle, the phase function can be determined by measuring the angular distribution of the scattered light in a goniometer.⁶ For this end, the diffusive sample is diluted such that single scattering events can be traced. However, in many cases, this straightforward approach is not possible. An important example is given by most biological tissues. Their phase function is determined by the spatial heterogeneities of the refractive index and thus depends among others on the cell contact surfaces that cannot be diluted without being destroyed. Typically, analyzed phase functions are characterized by the anisotropy value $g$ . If light propagation is modeled, often, a Henyey–Greenstein phase function (HGpf)⁷ with a specific anisotropy value is assumed. However, this procedure can give rise to errors up to 60% in the spatially resolved reflectance¹ for small source-detector separations. The anisotropy value does not describe the influence of the phase function on the reflectance properly. In order to resolve this issue, two parameters have been proposed to better describe the reflectance: Bevilacqua and Depeursinge⁸ proposed the principle of equivalence relations and its most prominent parameter $γ = (1 - g_{2}) / (1 - g_{1})$ with $g_{i}$ representing the $i$ ’th Legendre moment of the scattering phase function and thus $g_{1} = g$ the anisotropy value. This parameter is well suited to describe the reflectance in the spatial domain, however, not in the spatial frequency domain. To overcome this drawback, we recently introduced a new parameter

Eq. (1)

σ = \sum_{i = 0}^{\infty} {(- 0.5)}^{i} \frac{1 - g_{i + 2}}{1 - g_{1}},

for both the spatial and the spatial frequency domain.⁹ When the underlying phase function is not known in an experiment and a theory is fit to the data, a phase function parameter (e.g.,

g

,

γ

, or

σ

) can be adjusted in addition to

μ_{s}^{'}

and

μ_{a}

. The result can be best expressed in terms of

γ

or

σ

, since the latter parameters are little-dependent on the choice of phase function used in the theory as long as it covers the relevant range of

γ

or

σ

. The commonly used HGpf is not sufficient since it does not cover all the values

γ

or

σ

can reach. Chamot et al.¹⁰ stated a physiological range of

1.3 \leq γ \leq 2.4

corresponding to

0.74 \leq σ \leq 1.28

whereas the HGpf does not exceed

γ = 2

or

σ \approx 1.12

. Other models such as the Gegenbauer kernel phase function (GKpf) proposed by Reynolds and McCormick¹¹ with its parameters

α

and

g^{'}

are better suited in this regard. Kienle et al.² used a linear combination of two phase functions and fitted their ratios to retrieve

γ

. Turzhitsky et al.¹² applied a two parameter phase function derived from the Whittle-Matern correlation function to predict the subdiffusive reflectance. Kanick et al.⁴ provided a

γ

-map of a human hand via spatial frequency domain measurements. Very recently, Bravo et al.¹³ applied an inversion using a Monte Carlo (MC) lookup table and Naglic et al.¹⁴ proposed an inverse MC model based on

μ_{s}^{'}

,

μ_{a}

,

γ

, and the next order similarity parameter

δ

for describing the reflectance at very short source-detector separations. Despite these studies, the main reason not to regard the phase function in reflectance measurements remained the lack of an appropriate theory to describe the light propagation precisely and fast. The diffusion theory can be calculated very fast but does not take into account any phase function information. MC simulations can be applied to solve the radiative transfer equation (RTE); however, they need excessive computational time and are thus not applicable for real-time evaluations. Existing analytical solutions of the RTE¹⁵ can be used to fit a spatial reflectance curve. However, due to the additional phase function parameter, a tight incident beam compared to

1 / μ_{s}^{'}

and thus, high orders of the series expansion, coming with the need of calculations with quad double numbers, render the most difficult conditions for the analytical solution. As a result, fitting a single curve exceeds 30 min of single CPU core time. Thus, a faster method is required. Multivariate approaches are useful for this acceleration. In literature, artificial neural networks have been reported¹⁶^,¹⁷ to solve the inverse problem. Dam et al.¹⁸ reported minimal prediction errors using a principal component analysis (PCA), determining

μ_{s}^{'}

and

μ_{a}

in a small range of the possible optical parameters. Yaqin et al.¹⁹ improved their neural network using a PCA for preprocessing relative data. An interpolation from a lookup table with a large set of theoretical reflectance data is also possible but would require large capacity on a hard drive for three parameters to be varied. The reduction of this reflectance data can be performed by a PCA. This study presents a fast algorithm for the determination of

μ_{s}^{'}

,

μ_{a}

, and

σ

using an analytical solution of the RTE for the light propagation and a PCA for reducing the necessary data for a fast identification of the curve. This approach based on a three-parameter approximation improves the analysis of reflectance data for small source-detector separations and provides useful insights into the limits of the prediction and the modeling of the phase function, when realistic conditions are applied. When errors are given, they represent the median of the unsigned errors. For all the results in this study, the corresponding root mean square errors are, in average, 2.8 times (ranging from 1.4 to 8.4 for the individual examinations) larger than the median errors.

At first, the modeling and calculation of the light propagation according to the RTE are explained. In the next section, the principle of the evaluation is shown for only two parameters $μ_{s}^{'}$ and $μ_{a}$ , which is not only useful when the phase function is known but also much more demonstrative than the three parameter case. The method is then expanded to a three-parameter analysis. Section 3.2 is the main part of the report. It shows the capability for the determination of $μ_{a}$ , $μ_{s}^{'}$ , and the phase function parameter $σ$ in a very fast and robust way. The results from this new evaluation method are compared to those obtained by a fit.

2. Light Propagation Model

In this study, the light propagation is described by the nonpolarized RTE. We use the analytical derivation of the spatially resolved reflectance by Liemert and Kienle¹⁵ to exactly solve the $P_{N}$ -equations. We apply a Gaussian beam that illuminates the sample perpendicularly as shown in Fig. 1. The $P_{N}$ -equations are the result of the decomposition of the radiance into spherical harmonics. The decomposition into a high order $N$ is important in proximity to the source, where the radiance is a steep function in the angular space. However, both the high-order expansion and the small beam width come at the cost of increased computational time. Here, we mainly use $N = 11$ . Whenever we perform a fit that involves numerous calculations to find the minimum of the cost function, we reduce the order to $N = 7$ . Furthermore, one data set with $N = 17$ is calculated for comparison and to estimate the truncation error. The geometry of the sample is assumed to be the half space, and the boundary conditions are included exactly according to Fresnel’s law assuming a constant refractive index of $n_{M} = 1.40$ and $n = 1.00$ in- and outside the medium, respectively. In the analytical solution algorithm, also the phase function is projected onto spherical harmonics. These projection values are called “Legendre moments” $g_{i}$ . The GKpf¹¹ and its special case for $α = 0.5$ , the HGpf, are applied for most calculations in this study. In general, the corresponding Legendre moments are calculated numerically. For the HGpf, they can be calculated analytically using $g_{i} = g^{i}$ . Furthermore, we use a set of different types of phase functions described by Bodenschatz et al.⁹ For verification, we carry out MC simulations to solve the RTE using the same geometry except for the incident beam that is Dirac-shaped in the MC. This numerical solution is convolved with the same Gaussian beam profile assumed in the analytical calculations. After the convolution, the noise of the MC reflectance curves is negligible compared to the simulated measurement noise described below. This Gaussian beam decreases to $e^{- 2}$ at 0.25 mm corresponding to a full width half maximum of $\approx 0.3 mm$ . This tight beam profile enables a detailed investigation of the strong influence of the phase function close to the source but necessitates octuple-precision in the calculation of the analytical solution, which renders the calculation slow.

Fig. 1

Schematic of the experimental setup that corresponds to the calculated model.

In the study, we use sets of forward calculations of reflectance data that are described in Table 1. For the calculation of data set PCA1, used as data base for the 2D-PCA, we kept the phase function constant and calculated the spatially resolved reflectance for a grid of linearly spaced values of $μ_{s}^{'}$ and logarithmically spaced values of $μ_{a}$ in ranges indicated in the table. For the data base of the 3D-PCA (PCA2), for each point on the said $μ_{s}^{'} - μ_{a}$ -grid, the reflectance was calculated for 54 equally spaced values of $σ$ assuming a GKpf with $α = 1.8$ . The spacing of $σ$ was ensured using adequate values of the second GKpf parameter $g^{'}$ . The underlying optical coefficients of the test sets (TS A, TS B, TS C, TS D, and TS MC) were chosen randomly: for the given number of combinations, $μ_{s}^{'}$ was chosen as an equally distributed number in the range indicated in Table 1. The absorption value $μ_{a}$ is given by a random exponent in the presented range. The test sets differ in their choice of phase functions: in TS A, the constant phase function used for PCA 1 was used. In TS B, a GKpf was assumed with randomly chosen parameters $0 < α < 3$ and $- 0.5 < g^{'} < 1$ . If the corresponding $σ$ was not in the range $0.31 < σ < 1.29$ , the procedure was repeated. For TS C, the set of 1000 phase functions used in a former study of Bodenschatz et al.⁹ was used. For this purpose, for each combination of $μ_{s}^{'}$ and $μ_{a}$ , a random integer ranging from 1 to 1000 was chosen to indicate the phase function. In TS D, the phase function was chosen similar to TS B; however, $α$ was fixed to be 1.8. For TS MC, we assumed an HGpf and chose the anisotropy value $g$ in the range from 0.01 to 0.99. All data sets were calculated on a server machine with 16 cores except for the test set TS MC. The underlying MC simulations were calculated on the CPU of a desktop PC.

Table 1

Computational data overview. PN-curves represent analytical solutions of the RTE, MC simulation were used to solve the RTE numerically. x…y represents either a linear spacing (N steps) or the choice of a number from an equally distributed probability from x to y (random).

Name	Method	Curves	μs′·mm	μa·mm	Phase function
PCA 1	$P_{11}$	987	$0.5 \dots 5.1$ (47 steps)	$10^{- 2.8 \dots - 0.8}$ (21 steps)	HGpf, $g = 0.7$
PCA 2	$P_{11}$	12,960	$0.5 \dots 5.1$ (24 steps)	$10^{- 2.7 \dots - 0.9}$ (10 steps)	GKpf, $α = 1.8$ , $σ = 0.27 \dots 1.33$ (54 steps)
TS A	$P_{11}$	15,000	$0.63 \dots 4.2$ (random)	$10^{- 2.6 \dots - 1.0}$ (random)	HGpf, $g = 0.7$
TS B	$P_{11}$	8000	$0.63 \dots 4.2$ (random)	$10^{- 2.6 \dots - 1.0}$ (random)	GKpf: $0 < α < 3$ and $g^{'}$ random with $σ = 0.31 \dots 1.29$
TS C	$P_{11}$	15,000	$0.63 \dots 4.2$ (random)	$10^{- 2.6 \dots - 1.0}$ (random)	1000 phase functions, randomly chosen from Ref. 9
TS D	$P_{17}$	800	$0.63 \dots 4.2$ (random)	$10^{- 2.6 \dots - 1.0}$ (random)	GKpf: $α = 1.8$ , $g^{'}$ random with $σ = 0.31 \dots 1.29$
TS MC	MC	800	$0.63 \dots 4.2$ (random)	$10^{- 2.6 \dots - 1.0}$ (random)	HGpf, $g = 0.01 \dots 0.99$ (random)

We observe that it is very difficult to rely only on the form of the spatially resolved reflectance curve for the determination of the three optical properties $μ_{s}^{'}$ , $μ_{a}$ , and $σ$ . An absolute factor $K$ to scale the reflectance data $R$ has to be determined to compare the experimental and theoretical data. When this factor is known, we call the reflectance data “absolute,” otherwise, it is called “relative data.” In an experiment, the factor $K$ can be determined by calibrating the system using an optical phantom with known optical properties. Furthermore, curves called “ideal” have a factor $K$ that is exactly known. However, it is difficult to determine this factor precisely in an experiment. To evaluate the error in the determination of the optical properties caused by an error of $K$ , for curves called “realistic,” we multiply $K$ with a random factor between 0.95 and 1.05, so an error of up to 5% is assumed. Furthermore, for the “realistic” reflectance data $R_{l}$ (at points ${\tilde{ρ}}_{l} = l Δ x$ , $Δ x = 12.5 μ m$ representing the pixel separation in the sample plane), we assume normally distributed measurement noise with the standard deviation

Eq. (2)

Δ_{l} = \frac{R_{0}}{N_{e} \sqrt{2 π l}} \sqrt{N_{e} \frac{R_{l}}{R_{0}} + {(\sqrt{2} Δ_{RON})}^{2}},

that is added to the reflectance data

R_{l}

for each distance index

l

, which is at first a digital number given by the camera. Here,

N_{e}

represents the well capacity and

Δ_{RON}

the readout noise per pixel.

R_{0}

represents the digital number of the brightest pixel of the camera. The sum in the root consists of the shot noise and twice the readout noise, due to the necessity of taking a dark image. The dark current was neglected since it is much smaller than the readout noise for typical integration times less than 1s when using a cooled silicon-based camera. The fraction

R_{0} / N_{e}

scales the noise amplitude in electrons into the unit of

R_{l}

. It is divided by

\sqrt{2 π l}

since most camera-based spatially resolved reflectance setups integrate the signal over all pixels having the same distance to the center of the incident beam, thus, the number of pixels increases with the distance to the source. For the calculation, we set

N_{e} = 15,000

that represents 50% of the full well capacity and

Δ_{RON} = 1.2

. These values are given in the data sheet of the Zyla 5.5 (Andor, Belfast, Northern Ireland) to represent a typical scientific CMOS camera.

3. Inverse Problem Algorithms

3.1.

2D-PCA

In this section, the principle of the determination of $μ_{s}^{'}$ and $μ_{a}$ by use of a PCA (2D-PCA) is explained. A full tutorial on PCA can be found in Ref. 20. Briefly, PCA extracts the main data from a matrix in terms of orthogonal vectors. These vectors can be seen as the best-suited indexed family and thus concentrate the information of the matrix. In this study, only the main vectors are regarded, since for solving the inverse problem, the number of principal components should be at least equal to the number of parameters that are to be determined. To create the data matrix, the spatially resolved reflectance is binned such that it is given on a set of distances

Eq. (3)

ρ_{l} = β t_{l}^{2} + (ρ_{\max} - ρ_{\min} - β) t_{l} + ρ_{\min},

with

ρ_{\min} = 0.4 mm

,

ρ_{\max} = 6.0 mm

,

β = 5 mm

using

t_{l} = l / N_{L}

,

l = 1, \dots N_{L}

and the number of distances

N_{L} = 50

. This uncommon formalism of using not equally spaced positions increases the sensitivity to the phase function parameter, since there is a higher density of data points for smaller distances. In the result section, a comparison of this formalism with an equally spaced model is presented. Reflectance data for every combination

j

of the optical properties of the data set PCA 1 in Table 1 are calculated. To construct the matrix, on which PCA will be performed, we first define the matrix element

{\tilde{f}}_{j l}

with the reflectances

R

weighted with a power of the distance to counter the decrease of the reflectance signal

Eq. (4)

{\tilde{f}}_{j l} = ρ_{l}^{1.75} \cdot R [μ_{a, j}, μ_{s, j}^{'}, p (θ), ρ_{l}], l = 1, \dots, N_{L} and j = 1, \dots, 987 .

The 1.75th power of $ρ_{l}$ is a trade-off between the accuracies of the optical properties to be determined with a focus on $σ$ . A comparison of different weighting functions is shown below. For centralization and standardization of each row vector in this matrix, a mean absolute value $q_{0, j} = \sum_{l} {\tilde{f}}_{j l} / N_{L}$ called zeroth projection is extracted for the reflectance vector of each combination $j$ of optical properties. In order to get a variable that only depends on the shape of the reflectance curve, we extract the standard normal variate²¹ for each reflectance vector

Eq. (5)

f_{j l} = \frac{{\tilde{f}}_{j l} - q_{0, j}}{s_{j}},

with

s_{j}

representing the standard deviation of the weighted reflectance vector

{\tilde{f}}_{j l}

. Using this matrix

{(X)}_{j l} = f_{j l}

, we determine eigenvectors

{\vec{u}}_{m}

and the corresponding eigenvalues

λ_{m}

of its covariance matrix

X^{T} X

. The eigenvalues (together with their corresponding eigenvectors) are sorted by its values such that

λ_{1} \geq λ_{2} \geq λ_{3} \geq \dots

. The two eigenvectors

{\vec{u}}_{1}

and

{\vec{u}}_{2}

represent the relevant principal components. In the next step, the projection vectors

{\vec{q}}_{m} = X {\vec{u}}_{m}

are calculated for

m = 1,2

and normalized such that

⟨ {\vec{q}}_{m}, {\vec{q}}_{m} ⟩ = 1

with

m = 0,1, 2

. Please note that the vector

{\vec{q}}_{0}

was not obtained by the application of an eigenvector, but by calculating the mean absolute value. Nonetheless, it is very useful to align it with the projection values. Since every element of a projection vector

{\vec{q}}_{m}

represents a single combination of

μ_{a}

and

μ_{s}^{'}

, these vectors (

m = 0,1, 2

) are expanded onto a

μ_{a} - μ_{s}^{'}

-grid to become matrices

Q_{m}

with

m = 0,1, 2

. Using an analytical model, no noise impairs the data. This allows the smooth matrices

Q_{m}

to be interpolated using a finer grid to increase the resolution by a factor of two in both dimensions. An exemplary matrix

Q_{0}

is represented by the curved surface in Fig. 2(a). In summary, we have introduced the PCA algorithm to process the theoretical reflectance data. Now, a single reflectance data vector

\overset{ˇ}{\vec{R}}

given at positions

ρ_{l}

can be analyzed efficiently: processing

\overset{ˇ}{\vec{R}}

analogously, we get the mean absolute value

{\overset{ˇ}{Q}}_{0}

[horizontal plane in Fig. 2(a)] and the vector

\overset{ˇ}{\vec{f}}

. The projection values are given by the scalar product

{\overset{ˇ}{Q}}_{m} = ⟨ \overset{ˇ}{f}, {\vec{u}}_{m} ⟩

with

m = 1,2

. The optimal solution for

μ_{s}^{'}

and

μ_{a}

should minimize the differences of

{\overset{ˇ}{Q}}_{m}

and

Q_{m}

with

m = 0,1, 2

[see Figs. 2(b) and 2(c)]. We define the matrix of differences

ϕ = \sum_{m = 0}^{2} {(Q_{m} - {\overset{ˇ}{Q}}_{m})}^{2}

[displayed in Fig. 2(d)]. The discrete position

{\vec{x}}_{0}

in the grid yielding the minimal value of

ϕ

is analyzed. The true minimal position is assumed to be near

{\vec{x}}_{0}

and is determined by optimizing a quadric surface in the neighborhood of this value. This can be done very efficiently by calculating the Moore–Penrose pseudoinverse, which does not depend on the specific curve and thus only needs to be calculated once. The extremum of the quadric corresponds to the true minimal position that provides the optical properties. If the extremum is further away than the nearest neighbors on the grid,

{\vec{x}}_{0}

is assumed to represent the optical properties. The functions for solving the inverse problem by use of the PCA are implemented in Matlab2013b (The MathWorks Inc., Natick, Massachusetts).

Fig. 2

Principle of the determination of the two optical properties $μ_{a}$ and $μ_{s}^{'}$ . (a) The projection value of the curve to be analyzed (horizontal plane) is compared to that of curves of different optical properties. The square of the difference is taken [logarithmically displayed in (b) and for another projection value in (c)]. (d) These matrices are summed up. In the neighborhood of the lowest point, a quadric is used to find the optimal value to represent the determined optical properties (cross), compared to the optical properties of the forward calculation (circle).

3.2.

3D-PCA

For the determination of the three parameters $μ_{s}^{'}$ , $μ_{a}$ , and the phase function parameter $σ$ (3D-PCA), the procedures explained above are analogously applied. However, the transition of the vectors ${\vec{q}}_{m}$ to third-order tensors ( $μ_{a} - μ_{s}^{'} - σ$ -grid) instead of matrices ( $μ_{a} - μ_{s}^{'}$ -grid) renders a graphical demonstration difficult. The data set for finding the principal components consists of all combinations of optical properties described in PCA 1 in Table 1. Using a GKpf with $α = 1.8$ , the variation of $g^{'}$ corresponds to $σ$ ranging between 0.27 and 1.33. The projection values were interpolated on a grid with 72, 40, and 27 steps for $μ_{s}^{'}$ , $μ_{a}$ , and $σ$ , respectively. The resolution of the latter was reduced, since results with grids of higher resolution showed no relevant enhancement. Since three optical parameters are to be evaluated, the index $m$ in the calculation of $ϕ$ is extended to 3, the corresponding eigenvectors are presented in Fig. 3. Also, the result of the nonlinear spacing of $ρ$ can be found in the figure: since the eigenvectors are orthogonal, their oscillations have the highest frequency in the subdiffuse regime, where the influence of $σ$ is crucial.

Fig. 3

Four principal components $u_{i}$ , $i = 1,2, 3,4$ used in the 3D-PCA. The nonlinear spacing of $ρ$ can be seen in the increasing distance of the data points.

3.3.

Fit Routine

A nonlinear least squares optimization method²² is used to adjust the parameters $μ_{s}^{'}$ , $μ_{a}$ , $g^{'}$ as well as the scaling factor $K$ for relative data to minimize $\sum_{l} {[{\overset{ˇ}{R}}_{l} - K \cdot R (μ_{s}^{'}, μ_{a}, g^{'}, ρ_{l})]}^{2} / Δ_{l}$ with evenly spaced distances from the source $ρ_{l}$ . The high precision and corresponding slow computation render the use of the $P_{11}$ -solution for curve fitting difficult. We therefore make use of the $P_{7}$ -solution that allows us to perform computations of a large and statistically meaningful data set. A GKpf with $α = 1.8$ is assumed, motivated by the large range of $σ$ possible by varying $g^{'}$ .

4. Results

4.1.

2D-PCA

In this section, a large number of reflectance curves are analyzed. This number is given in the corresponding rows in Table 1. At first, test set A (TS A) is analyzed by the 2D-PCA. The analysis of a single curve runs on one CPU core in less than 1 ms after initial loading of the prepared data. The determined optical properties $μ_{s}^{'}$ and $μ_{a}$ for ideal and realistic curves (see Sec. 2) are shown in Figs. 4(a) and 4(b), respectively. It can be found that the predicted optical properties of the ideal curves match the expected values almost perfectly. The deviations of the obtained optical properties for the realistic curves give an impression on how distorted the information is. Please note that in this step the phase function is assumed to be known and thus used in the forward calculations of PCA 1. To statistically demonstrate the quality of the optical property determination for the ideal curves, boxplots are shown in Fig. 5 featuring the 25% and the 75% quartiles as well as the median values of the relative error of each optical property. The whiskers usually displayed to represent the 5% and the 95% quantiles are not shown since they are mostly two to six times the distance to the median compared to that of the corresponding quartiles and thus confirm the conclusions drawn from the quartiles. Only for the small errors, this factor is exceeded for six whiskers. In the visualization using the boxplots, the perfect match for the ideal curves of TS A can be found too. The sensitivity of the 2D-PCA algorithm to an unknown phase function is shown in Fig. 6, where the 2D-PCA of TS B is presented. This additional degree of freedom makes it more difficult to determine $μ_{s}^{'}$ and $μ_{a}$ . The errors in $μ_{s}^{'}$ and $μ_{a}$ caused by a phase function difference exceed those caused by the realistic noise assumptions and distortion of $K$ by far. Thus, for the two parameter determination of $μ_{s}^{'}$ and $μ_{a}$ , we show the analysis of ideal curves by the 2D-PCA in Figs. 5 and 6 and do not show the similar results of the realistic curves. Since $σ$ was not determined in that case, in the right column, no corresponding box is shown for the 2D-PCA.

Fig. 4

Optical properties (a) $μ_{s}^{'}$ and (b) $μ_{a}$ of TS A determined by the use of 2D-PCA. Each point represents the solution of one forward curve of TS A.

Fig. 5

Boxplots quantifying errors in the determination of the optical properties using 2D- and 3D-PCA. The 25%, 50% (median), and 75% quantiles are shown. Three boxes are drawn together that represent the analysis of the ideal curves by 3D-PCA (top), the realistic curves by 3D-PCA (center), and the ideal curves by 2D-PCA (bottom) for the parameters (a) $μ_{s}^{'}$ , (b) $μ_{a}$ , and (c) $σ$ . The 5% and the 95% quantiles (not shown) are two to four times the distance to the median than the corresponding quartile (one exception).

Fig. 6

Optical properties of TS B determined by the use of 2D-PCA. The results for the realistic curves are not shown because of their similarity to the results of the ideal curves.

For the other test sets, similar results can be found in Fig. 5. Analyzing TS C, different kinds of phase functions produce large deviations in $μ_{s}^{'}$ and $μ_{a}$ . Interestingly, in the analysis of TS MC, the range of deviations is smaller compared to the other test sets with varying phase functions. This is probably due to the fact that only HGpf is used for the MC test set as in the forward calculation of PCA 1. It underlines the limited variability of HGpf even if $g$ is varied in a physically reasonable range.

4.2.

3D-PCA

Next, the 3D-PCA is applied on the test sets. Because of the additional dimension, the computational effort is increased by about one order of magnitude. The determination of the optical properties for a single curve takes around 10 ms of CPU time on a single core. When TS B is analyzed by the 3D-PCA [see Figs. 7(a)–7(c)], the prediction of $μ_{s}^{'}$ and $μ_{a}$ is improved for both ideal and realistic curves compared to the 2D-PCA. The results for the ideal curves represent the errors of the 3D-PCA algorithm. Those of the realistic curves show the influence of the noise assumption. The study of Fig. 7 reveals some horizontal and vertical clustering of data points. They correspond to solutions that could not be optimized using the quadric fitting (see Sec. 3.1) and thus lay on the grid points of the forward calculations. The determined values of the phase function parameter $σ$ are in good agreement with the expected values for ideal curves. For realistic curves, however, its determination error increases rapidly. This is not due to the added noise that has only a small effect on the determination of the optical properties, but mainly due to the imprecisely known absolute value $K$ . The increase of the errors of the realistic case compared to the ideal one is visualized more clearly in Fig. 5. For realistic curves, 50% of the predicted values of $σ$ are less than 0.027 from the expected value corresponding to 3% of the full range of $σ$ ( $0.31 \leq σ \leq 1.29$ ) and 5% of the physiologically relevant values ( $0.74 \leq σ \leq 1.28$ ). The latter values are obtained by plotting $σ$ over $γ$ for the 1000 phase functions used in TS C and extracting the “physiological range” of $1.3 \leq γ \leq 2.4$ .¹⁰ An important and most difficult challenge is the determination of the optical properties in TS C by 3D-PCA [see Figs. 7(d)–7(f)], since the underlying phase function is not known at all. In the corresponding boxplots in Fig. 5(a), a small systematic underestimation of about 1% can be found for $μ_{s}^{'}$ . The uncertainties of $μ_{a}$ are comparable to those of TS B. In the phase function parameter determination, the large variation of the phase functions used for TS C results in the largest median errors of $σ$ represented by the boxes. For ideal curves, the influence of this variation can be seen most clearly. The median is shifted by about 0.015 and the quartiles extend to $\sim 0.03$ and are thus most distant to the median compared to any other test set used in the study. For realistic curves, the determination is only slightly worse compared to TS B. The analysis of TS D ( $P_{17}$ ) illustrates the effect of the limited order $N = 11$ in the calculation of PCA 2. Compared to TS B, no significant difference in the determination of $μ_{s}^{'}$ and $μ_{a}$ can be found. In the determination of $σ$ , a small drift of the median of less than 0.01 can be found for the analysis of realistic curves. Since TS D uses the same type of phase functions as PCA 2, this drift can only be explained with their different orders $N$ . Since both $σ$ and the calculation order $N$ have an increased influence on the reflectance at small distances from the source, the more precise description of the reflectance from $P_{17}$ is falsely interpreted as a change in $σ$ . However, the observed deviation in $σ$ is smaller than the typical errors even for ideal curves. That means, a calculation order of $N = 11$ is sufficient. Finally, the 3D-PCA is applied on reflectance curves from MC simulations. The results for $μ_{s}^{'}$ and $μ_{a}$ are comparable to the results from the TS B analysis, only the deviation in $σ$ is larger. This is probably due to the difference in the phase function type in TS MC (GKpf with $α = 0.5$ ) and in PCA 2 (GKpf with $α = 1.8$ ).

Fig. 7

Optical properties of (a–c) TS B and (d–f) TS C determined by the use of 3D-PCA.

4.3.

Configurations

The parameters that are used in the evaluations are not unique. In the following, we present results of the analysis when parameters are changed. The analysis of TS C using the 2D-PCA and the 3D-PCA is performed as a metric. Although the 2D-PCA could easily produce better results by choosing other parameters, it is the aim of this study to present the 3D-PCA and for simplicity, we choose to keep the parameters identical for both methods. In Fig. 8, the results of $μ_{a}$ and $σ$ are shown for different variations of the underlying parameters. At first, we demonstrate the inability of the determination of $σ$ using a linear spacing of $ρ$ compared to the nonlinear binning presented in Eq. (3). Using the linear representation of $ρ$ , the determination of $σ$ is strongly biased and shows larger errors than in the nonlinear case. Furthermore, the accuracy of the determination of $μ_{a}$ using the 3D-PCA is significantly improved using the nonlinear representation, whereas the results of $μ_{a}$ using the 2D-PCA are comparable. The results of $μ_{s}^{'}$ are less than 10% for both spacings and all examined weightings and allow similar conclusions as the discussed results of $μ_{a}$ and are thus not shown. The best weighting (reference) was chosen to be $ρ^{1.75}$ . It allows a determination of $σ$ with a high accuracy without increasing the error of the 2D-PCA too much. Please remark that the chosen parameters are to be optimized for the 3D-PCA; however, very large errors in the 2D-PCA indicate a nonrobust analysis. The commonly used linear spacing and weighting of $\log [R (ρ)]$ shows the smallest errors for $μ_{a}$ combined with good results for $μ_{s}^{'}$ in the 2D-PCA. It is, however, not suited for the 3D-PCA as it produces large errors in $μ_{a}$ and $σ$ . In this study, the number of used projection values is one larger than the number of optical properties to be determined. Using more projection values, the determination of $σ$ gets slightly more accurate at the cost of a strong deterioration of $μ_{a}$ and additionally in the 2D-PCA of $μ_{s}^{'}$ . To use one projection value less works equally precise in the 2D-PCA but deteriorates the results in the 3D-PCA for all three parameters. The number of data point was varied between 20 and 100. This modification, however, is not shown since it did not produce any difference. Lower numbers of data points are interesting for fiber probes that are not examined in this study.

Fig. 8

Comparison of variations in the analysis: the result of 2D-PCA is shown in terms of the error of prediction of the absorption coefficient of TS C. The 3D-PCA is represented by the errors in the determination of $μ_{a}$ and $σ$ of TS C. The configuration “reference” is used for all other results presented in this study and is chosen as a trade-off with focus on the determination of $σ$ .

4.4.

Least Square Fit

At last, we compare the presented methods to a fit algorithm. 1000 curves out of TS C were fit by the analytical $P_{7}$ -solution via variation of $μ_{s}^{'}$ , $μ_{a}$ , $K$ , and $g^{'}$ . A GKpf was assumed using $α = 1.8$ . Three different modalities were examined: in the first (a), the start parameters for the fit were $μ_{s}^{'} = 2 {mm}^{- 1}$ , $μ_{a} = 0.01 {mm}^{- 1}$ , $K = 1$ , $g^{'} = 0.5$ , in the second (b), the start parameters were given by the result of the 3D-PCA. In both (a) and (b), the theory was fit to ideal curves, whereas in the third (c), realistic curves were fit with the 3D-PCA start parameters as used in (b). The results of these fits are presented as boxplots in Fig. 9. The 5%- and 95%-whiskers are not shown; however, this time because they are too large and exceed the axes’ limits by far.

Fig. 9

Boxplots quantifying errors of the fit results of TS C in different configurations. (a) Fits of ideal curves with constant start parameters and (b) with start parameters given by the result of the 3D-PCA were performed, additionally, (c) realistic curves were fit using start parameters given by the 3D-PCA.

It can clearly be seen that the fit algorithm enhances the prediction results compared to the 3D-PCA, most noticeably for the parameter $μ_{a}$ , where the limits of the graph were copied from Fig. 5. Interestingly, modality (c) shows no deterioration even though realistic curves were analyzed instead of ideal ones. In the relative fit, the scattered absolute intensity can be compensated for by variation of $K$ , so only the effect of noise remains. Summarized, the prediction of $μ_{s}^{'}$ and $μ_{a}$ is enhanced, the prediction of $σ$ is slightly worse compared to the 3D-PCA of TS C with ideal curves and gives better results for the realistic ones. The main difference and the reason for the enhancement of $μ_{s}^{'}$ and $μ_{a}$ is the weighting in the fit that uses the noise applied to achieve realistic curves. That measure shifts the focus from $σ$ to $μ_{a}$ . In the 3D-PCA, the weighting is done by choosing nonequidistant values of $ρ$ , which is justified but arbitrary and by multiplication of the reflectance with $ρ^{2}$ in Eq. (4). Concerning $σ$ , the systematic deviation in the fit is partly due to the choice of phase function in TS C (in 3D-PCA there is also a drift, but smaller) and partly due to the lower calculation order of $N = 7$ in the fit. The choice of a lower order is unfortunate but necessary. The computational effort for the fit is immense, even using $N = 7$ . An average total CPU time of 34.0, 14.7, and 14.2 min was needed per curve using the modalities (a) to (c), respectively. Please note that the fit algorithm can severely miss the expected optical properties, these values cannot be seen in the quartiles. Hence, only the best 50% of the results of each algorithm are compared.

In order to evaluate the validity of geometrical approximations, MC simulations were performed regarding an obliquely incident beam tilted by 5 deg and a numerical aperture angle of 5 deg of the detection. These variations were included both separately and at the same time. The optical properties were $μ_{s}^{'} = 1 {mm}^{- 1}$ , $μ_{a} = 0.01 {mm}^{- 1}$ , $n = 1.40$ , and $g = 0.7$ using an HGpf. The resulting curves were fit by the analytical $P_{11}$ solution assuming perpendicular incidence and full numerical aperture of the detection system. An error of 2%, 10%, and 0.01 was not exceeded for any of these cases in $μ_{s}^{'}$ , $μ_{a}$ , and $σ$ , respectively.

5. Conclusion

We presented a fast and robust method based on PCA to determine the reduced scattering coefficient $μ_{s}^{'}$ , the absorption coefficient $μ_{a}$ , and the recently introduced phase function parameter $σ$ from spatially resolved reflectance data. When these parameters and the underlying type of phase functions are not known and the reflectance data is not impaired with noise, we found the median errors of these parameters to be $Δ μ_{s}^{'} = 1.4 %$ , $Δ μ_{a} = 5 %$ , and $Δ σ = 0.025$ . The assumption of realistic measurement uncertainties increases these errors to $Δ μ_{s}^{'} = 2.5 %$ , $Δ μ_{a} = 12 %$ , and $Δ σ = 0.04$ . After preparation of the principal components, the calculation time is about 10 ms for the analysis of a single curve using MATLAB^® on a single core on a PC built in 2011. Absolute reflectance data are needed for successful determination of the optical properties, the influence of errors up to 5% in this absolute value was examined. If the phase function is known, the determination of $μ_{s}^{'}$ and $μ_{a}$ is highly accurate, takes only 1 ms per curve, and can even be achieved using relative data. Compared to a fit algorithm, the computational effort is decreased by orders of magnitude; however, the speedup comes at the price of a lower accuracy. If very high accuracy is needed, a further application of our method is the use of the 3D-PCA to find the start values for the fit to reduce its computational time by a factor of 2 to 4 while maintaining the high prediction quality of the fit routine and even increasing its robustness. If the phase function is kept constant in the determination of $μ_{s}^{'}$ and $μ_{a}$ from a set of ideal reflectance curves calculated for multiple phase functions, median errors of 5% in $μ_{s}^{'}$ and 22% in $μ_{a}$ occur. Hence, even for cases, where information about the phase function is not of interest or relevance, the presented method improves the determination of $μ_{s}^{'}$ and $μ_{a}$ by giving the possibility to adjust the phase function. Neither the presented algorithm nor the fit routine was able to determine $σ$ better than $Δ σ = 0.04$ . These uncertainties in $σ$ relate to the limits of describing phase function influence on spatially resolved reflectance by a single parameter. The calculation of the data set for the analysis (PCA 2) can thus be accelerated by assuming only 10 to 15 different values for $σ$ without reducing the determination accuracy. Please note that $γ$ could also be used as a reproducible phase function parameter in the spatial domain. However, in view of the more universal applicability of $σ$ also to the spatial frequency domain, $σ$ allows for enhanced comparability.

Future work includes the analysis of the experimental data and the implementation of the PCA into the automated measurement software to allow for real-time extraction of optical properties from spatially resolved reflectance measurements.

Disclosures

No conflict of interest, financial or otherwise, are declared by the authors.

Acknowledgments

The authors thank the reviewers for their comprehensive comments that helped to improve the final results of this study. P. Krauter acknowledges the support by Evangelisches Studienwerk Villigst e.V.

References

1.

J. R. Mourant et al., “Influence of the scattering phase function on light transport measurements in turbid media performed with small source–detector separations,” Opt. Lett., 21 (7), 546 –548 (1996). http://dx.doi.org/10.1364/OL.21.000546 OPLEDP 0146-9592 Google Scholar

2.

A. Kienle, F. K. Forster and R. Hibst, “Influence of the phase function on determination of the optical properties of biological tissue by spatially resolved reflectance,” Opt. Lett., 26 (20), 1571 –1573 (2001). http://dx.doi.org/10.1364/OL.26.001571 OPLEDP 0146-9592 Google Scholar

3.

N. Bodenschatz et al., “Detecting structural information of scatterers using spatial frequency domain imaging,” J. Biomed. Opt., 20 (11), 116006 (2015). http://dx.doi.org/10.1117/1.JBO.20.11.116006 JBOPFO 1083-3668 Google Scholar

4.

S. C. Kanick et al., “Sub-diffusive scattering parameter maps recovered using wide-field high-frequency structured light imaging,” Biomed. Opt. Express, 5 (10), 3376 –3390 (2014). http://dx.doi.org/10.1364/BOE.5.003376 BOEICL 2156-7085 Google Scholar

5.

P. Naglič et al., “Estimation of optical properties by spatially resolved reflectance spectroscopy in the subdiffusive regime,” J. Biomed. Opt., 21 (9), 095003 (2016). http://dx.doi.org/10.1117/1.JBO.21.9.095003 JBOPFO 1083-3668 Google Scholar

6.

F. Foschum and A. Kienle, “Optimized goniometer for determination of the scattering phase function of suspended particles: simulations and measurements,” J. Biomed. Opt., 18 (8), 085002 (2013). http://dx.doi.org/10.1117/1.JBO.18.8.085002 JBOPFO 1083-3668 Google Scholar

7.

L. Henyey and J. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J., 93 70 –83 (1941). http://dx.doi.org/10.1086/144246 Google Scholar

8.

F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source–detector separations close to one transport mean free path,” J. Opt. Soc. Am. A, 16 (12), 2935 –2945 (1999). http://dx.doi.org/10.1364/JOSAA.16.002935 JOAOD6 0740-3232 Google Scholar

9.

N. Bodenschatz et al., “Quantifying phase function influence in subdiffusively backscattered light,” J. Biomed. Opt., 21 (3), 035002 (2016). http://dx.doi.org/10.1117/1.JBO.21.3.035002 JBOPFO 1083-3668 Google Scholar

10.

S. Chamot et al., “Physical interpretation of the phase function related parameter γ studied with a fractal distribution of spherical scatterers,” Opt. Express, 18 (23), 23664 –23675 (2010). http://dx.doi.org/10.1364/OE.18.023664 OPEXFF 1094-4087 Google Scholar

11.

L. Reynolds and N. McCormick, “Approximate two-parameter phase function for light scattering,” J. Opt. Soc. Am., 70 (10), 1206 –1212 (1980). http://dx.doi.org/10.1364/JOSA.70.001206 JSDKD3 Google Scholar

12.

V. Turzhitsky et al., “A predictive model of backscattering at subdiffusion length scales,” Biomed. Opt. Express, 1 (3), 1034 –1046 (2010). http://dx.doi.org/10.1364/BOE.1.001034 BOEICL 2156-7085 Google Scholar

13.

J. J. Bravo et al., “Sub-diffuse optical biomarkers characterize localized microstructure and function of cortex and malignant tumor,” Opt. Lett., 41 (4), 781 –784 (2016). http://dx.doi.org/10.1364/OL.41.000781 OPLEDP 0146-9592 Google Scholar

14.

P. Naglič et al., “Adopting higher-order similarity relations for improved estimation of optical properties from subdiffusive reflectance,” Opt. Lett., 42 (7), 1357 –1360 (2017). http://dx.doi.org/10.1364/OL.42.001357 OPLEDP 0146-9592 Google Scholar

15.

A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep., 3 2018 (2013). http://dx.doi.org/10.1038/srep02018 SRCEC3 2045-2322 Google Scholar

16.

T. Farrell, B. Wilson and M. Patterson, “The use of a neural network to determine tissue optical properties from spatially resolved diffuse reflectance measurements,” Phys. Med. Biol., 37 (12), 2281 –2286 (1992). http://dx.doi.org/10.1088/0031-9155/37/12/009 PHMBA7 0031-9155 Google Scholar

17.

M. Jäger, F. Foschum and A. Kienle, “Application of multiple artificial neural networks for the determination of the optical properties of turbid media,” J. Biomed. Opt., 18 (5), 057005 (2013). http://dx.doi.org/10.1117/1.JBO.18.5.057005 JBOPFO 1083-3668 Google Scholar

18.

J. S. Dam et al., “Determination of tissue optical properties from diffuse reflectance profiles by multivariate calibration,” Appl. Opt., 37 (4), 772 –778 (1998). http://dx.doi.org/10.1364/AO.37.000772 APOPAI 0003-6935 Google Scholar

19.

C. Yaqin et al., “Determination of tissue optical properties from spatially resolved relative diffuse reflectance by PCA-NN,” in Neural Networks and Signal Processing, 2003. Proc. of the 2003 Int. Conf. on, 369 –372 (2003). http://dx.doi.org/10.1109/ICNNSP.2003.1279286 Google Scholar

20.

L. Smith, “A tutorial on principal component analysis,” (2014) http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf September ). 2014). Google Scholar

21.

R. Barnes, M. S. Dhanoa and S. J. Lister, “Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra,” Appl. Spectrosc., 43 (5), 772 –777 (1989). http://dx.doi.org/10.1366/0003702894202201 APSPA4 0003-7028 Google Scholar

22.

S. Agarwal et al., “Ceres solver,” (2015) http://ceres-solver.org April ). 2015). Google Scholar

Biographies for the authors are not available.

Citation Download Citation

Philipp Krauter, Dominik Reitzle, Simeon Geiger, and Alwin Kienle "Determination of three optical properties from subdiffusive spatially resolved reflectance calculations," Journal of Biomedical Optics 22(7), 075003 (21 July 2017). https://doi.org/10.1117/1.JBO.22.7.075003

Received: 23 March 2017; Accepted: 5 July 2017; Published: 21 July 2017

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY