## Abstract

This paper is the second of two dealing with light diffusion in a turbid cylinder. The diffusion equation was solved for an N-layered finite cylinder. Solutions are given in the steady-state, frequency, and time domains for a point beam incident at an arbitrary position of the first layer and for a circular flat beam incident at the middle of the cylinder top. For special cases the solutions were compared to other solutions of the diffusion equation showing excellent agreement. In addition, the derived solutions were validated by comparison with Monte Carlo simulations. In the time domain we also derived a fast solution (≈ 10ms) for the case of equal reduced scattering coefficients and refractive indices in all layers.

©2010 Optical Society of America

## 1. Introduction

In the accompanying paper we described the light diffusion in a homogeneous cylindrical turbid medium [1]. In this study we present solutions of the diffusion equation for cylindrical turbid media having an arbitrary number of layers.

The light propagation in layered turbid media has been investigated for therapeutical or diagnostic applications because many parts of the body exhibit a layered structure, e.g. the arm, the leg, or the head. The transport theory and an approximation of it, the diffusion theory, have been normally applied to obtain quantitative data [2]. The advantage of the diffusion equation is that it can be solved analytically for layered geometries yielding fast and accurate solutions compared to the normally applied numerical solutions.

For laterally infinitely extended turbid media Dayan *et al*. presented approximative solutions of the two-layered diffusion equation [3]. We derived exact solutions of the two-layered diffusion equation in the steady-state, frequency, and time domains using the Fourier transform formalism [4, 5, 6]. This approach was, recently, extended by us and others for an infinite number of layers [8, 7, 9, 10], whereas the diffusion equation for a finite number of layers was investigated also by several other groups [11, 12, 13, 14, 15, 16].

A solution of the layered diffusion equation for a cylindrical geometry has been presented by Martelli *et al*. for the case of two and three layers [17, 18, 19]. They used the eigenfunction method together with the separation of variables to derive solutions in the time domain.

In this study we present solutions of the diffusion equation for an N-layered cylinder. The solutions are given in the steady-state, frequency, and time domains. For equal reduced scattering coefficients and refractive indices in all layers we also derived a simplified solution in the time domain. For special cases, for which solutions of the diffusion equation are described in literature, we compared these solutions to those derived in this study. In addition, we validated the derived solutions with Monte Carlo simulations.

Besides in vivo applications the derived solutions can be used especially for experiments on layered phantoms, because they often have a cylindrical geometry. Thus, the phantoms’ volume does not have to be so large that the assumption of a semi-infinite geometry is valid. In addition, the obtained solutions can be easily transferred to solutions of the heat conduction equation, where multi-layered cylindrical geometries are often encountered.

## 2. Theory

#### 2.1. Diffusion Theory

Solutions of the diffusion theory for an N-layered cylinder are, first, derived in the frequency domain, which include for the special case of zero modulation frequency the solutions in the steady-state domain. Then, solutions in the time domain are presented. The solutions are derived for a point beam incident onto an arbitrary location of the first layer. Thus, the beam can be directed onto the cylinder top as well as onto the cylinder barrel of the first layer. Additionally, we give the solutions for a circular flat beam incident onto the middle of the cylinder top.

Fig. 1 shows a scheme of the considered N-layered cylinder. The thickness, the refractive index, the reduced scattering and absorption coefficients of layer k are denoted by *l _{k}*,

*n*,

_{k}*μ*′

*, and*

_{sk}*μ*, respectively. The radius of the cylinder is indicated with

_{ak}*a*. As usual, it is assumed that the incident beam can be represented by an isotropic source at a distance of 1/(

*μ*′

_{s1}+

*μ*

_{a1}) from the location of incidence at the boundary of the cylinder. The point source is located in the first layer of the cylinder.

As in the accompanying paper the extrapolated boundary condition is used for the top, bottom, and barrel boundaries [1]. The extrapolation length is calculated with

where *R*
_{eff,k} is the angle-averaged probability for reflection at the boundary between layer k and the surrounding medium and *D _{k}* = 1/(3

*μ*′

_{sk}) is the diffusion coefficient of layer k. For the boundary at the top and the bottom of the cylinder k equals 1 and N, respectively. For the boundary at the cylinder barrel k is assumed to be the layer at which the remitted or transmitted light is calculated, because it is more involved to solve the diffusion equation for N layers using a different boundary lengths at the cylinder barrel for each layer.

### 2.1.1. Solution in the frequency domain

The diffusion equation in the frequency domain is given by

where Φ, *c*, and *ω* denote the fluence rate, the speed of light, and the angular frequency of the intensity modulated light, respectively. The source function for a point source in cylindrical coordinates is $S(\overrightarrow{r},\omega )=\frac{1}{\rho}\delta \left(\rho -{\rho}_{0}\right)\delta \left(\varphi -{\varphi}_{0}\right)\delta \left(z-{z}_{0}\right)$. Eq. 2 can be rewritten by using cylindrical coordinates as follows:

In the accompanying paper [1] we presented an integral transform, which can be used to simplify Eq. 3 to an ordinary differential equation. Due to the extrapolated boundary condition this transform is used in the following form

where *s _{n}* are determined from the following equation

*a*′ = *a* + *z _{bk}*, and

*J*is the first kind Bessel function of order m. By applying Eq. 4 to Eq. 3 the resulting ordinary differential equation for Φ = Φ(

_{m}*s*,

_{n}*ν*,

*m*,

*z*,

*ω*) is

The solution of Eq. 6 is obtained by solving

where $\alpha =\sqrt{{\mu}_{a}/D+{s}_{n}^{2}+\mathrm{i\omega}/\mathrm{Dc}}$, and by using

For the first layer the solution of Eq. 7 is composed of the homogeneous and particular part, for the other layers only by the homogeneous part. For deriving the particular solution of the first layer *G*
_{1}
^{(p)}(*s _{n}*,

*z*,

*ω*) we use the Fourier transform to find the one-dimensional Green’s function for an unbounded region in

*z*-direction. By applying

to Eq. 7 we get the Green’s function in the Fourier-space as

The inverse Fourier transform

is used to obtain the particular solution of Eq. 7 for *k* = 1

Thus, the complete solution of Eq. 7 is

$${G}_{k}({s}_{n},z,\omega )={A}_{k}{e}^{{\alpha}_{k}z}+{B}_{k}{e}^{-{\alpha}_{k}z},\phantom{\rule{5em}{0ex}}k=2,\dots ,N.$$

The constant coefficients *A _{k}* und

*B*are determined by using the following boundary conditions in z-direction

_{k}$${n}_{k+1}^{2}{G}_{k}\left({s}_{n},z={L}_{k},\omega \right)={n}_{k}^{2}{G}_{k+1}\left({s}_{n},z={L}_{k},\omega \right)$$

$${D}_{k}\partial {G}_{k}\left({s}_{n},z={L}_{k},\omega \right)/\partial z={D}_{k+1}\partial {G}_{k+1}\left({s}_{n},z={L}_{k},\omega \right)/\partial z$$

$${G}_{N}\left({s}_{\mathrm{nN}},z={L}_{N}+{z}_{\mathrm{bN}},\omega \right)=0,$$

where ${L}_{k}={\displaystyle \sum _{i=1}^{k}}{l}_{i},\phantom{\rule{1em}{0ex}}k=1,\dots ,N$. The solution for the fluence rate in the first layer (0 ≤ *z* < *l*
_{1}) is given by

$$\phantom{\rule{7em}{0ex}}\times \frac{{D}_{1}{\alpha}_{1}{n}_{1}^{2}{\beta}_{3}-{D}_{2}{\alpha}_{2}{n}_{2}^{2}{\gamma}_{3}}{{D}_{1}{\alpha}_{1}{n}_{1}^{2}{\beta}_{3}\mathrm{cosh}\left[{\alpha}_{1}\left({l}_{1}+{z}_{b1}\right)\right]+{D}_{2}{\alpha}_{2}{n}_{2}^{2}{\gamma}_{3}\mathrm{sinh}\left[{\alpha}_{1}\left({l}_{1}+{z}_{b1}\right)\right]}$$

and for the fluence rate in the N^{th} layer (*L*
_{N-1} ≤ *z* < *L _{N}*) by

In general, the quantities *β*
_{3} and *γ*
_{3} are obtained by recurrence relations. The start values are

$$\phantom{\rule{2em}{0ex}}+{D}_{N}{\alpha}_{N}{n}_{N}^{2}\mathrm{sinh}\left({\alpha}_{N-1}{l}_{N-1}\right)\mathrm{cosh}\left[{\alpha}_{1}\left({l}_{N}+{z}_{b2}\right)\right]$$

$${\gamma}_{N}={D}_{N-1}{\alpha}_{N-1}{n}_{N-1}^{2}\mathrm{sinh}\left({\alpha}_{N-1}{l}_{N-1}\right)\mathrm{sinh}\left[{\alpha}_{N}\left({l}_{N}+{z}_{b2}\right)\right]$$

$$\phantom{\rule{2em}{0ex}}+{D}_{N}{\alpha}_{N}{n}_{N}^{2}\mathrm{cosh}\left({\alpha}_{N-1}{l}_{N-1}\right)\mathrm{cosh}\left[{\alpha}_{N}\left({l}_{N}+{z}_{b2}\right)\right]$$

and the recurrence relations are given by

$${\gamma}_{k-1}={D}_{k-2}{\alpha}_{k-2}{n}_{k-2}^{2}\mathrm{sinh}\left({\alpha}_{k-2}{l}_{k-2}\right){\beta}_{k}+{D}_{k-1}{\alpha}_{k-1}{n}_{k-1}^{2}\mathrm{cosh}\left({\alpha}_{k-2}{l}_{k-2}\right){\gamma}_{k}.$$

For a two-layered cylinder (N=2) Eqs. 17 and 18 are not needed, instead *β*
_{3} = sinh[*α*
_{2}(*l*
_{2}+*z*
_{b2})] und *γ*
_{3} = cosh[*α*
_{2}(*l*
_{2}+*z*
_{b2})] have to be applied. For a three-layered cylinder (N = 3) only Eq. 17 is required. We note that we derived the above presented solutions already for the case of a laterally infinitely extended medium [9], but now they are presented in a different way, so that no approximations are needed anymore for avoiding numerical errors as was necessary for the previous solutions. Thus, we recommend to use Eqs. 15–18 also for the laterally infinitely extended medium. The inverse relation of Eq. 4 is given by [1]

Finally, the solution of Eq. 3 in the frequency domain for layer *k* is obtained using Eq. 8 and the following recurrence relation

as

For a point source that is incident onto the center of the cylinder top Eq. 21 is simplified to

For the derivation of the solutions of the N-layered diffusion equation for a circular flat beam with radius *ρ _{w}* incident onto the middle of the cylinder top we refer to the accompanying paper [1] where we gave the corresponding solution for a homogeneous cylinder. Accordingly, the solution of an N-layered cylinder is

We note that the solutions in the steady-state domain can be easily obtained from the above presented solutions in the frequency domain by setting *ω*= 0.

### 2.1.2. Solution in the time domain

In order to obtain the solutions of the diffusion equation in the time domain

we use the inverse Fourier transform

Thus, we can apply the solutions derived in the frequency domain to obtain the solutions in the time domain. For example, for a point source incident at an arbitrary position onto the first layer of the N-layered cylinder we obtain by inserting Eq. 21 into Eq. 25

The inverse Fourier integral is calculated by using the inverse discrete Fourier transform. For equal reduced scattering coefficients and refractive indices in all N layers faster solutions can be derived in the time domain. For this case, the complex numbers *α _{k}* in the frequency domain have the form

The definition of the Laplace transform with the complex frequency *s* = *σ*+*iω* is

where ROC is the region of convergence. By choosing *σ*=*μ _{ak}c*+

*Dcs*

^{2}

*one obtains*

_{n}Now it can be seen that Eq. 29 describes the Laplace transform of the function *x*(*t*)*e*
^{-Dcsn2t}. The region of convergence for this case is Re{*s*} > -*μ _{ak}c*. The imaginary axis is within the region of convergence, so that we can calculate the inverse transform of

*X*(

*μ*+

_{a}c*iω*) by means of the Fourier transform, in particular the FFT algorithm. For this special case we obtain the solution as

Finally, we give the results for a point source incident onto the center of the cylinder top

We note that the solutions for equal reduced scattering coefficients and refractive indices in each layer can be quickly evaluated (≈ 10ms) using the above equations.

The derived solutions of the diffusion equation give the fluence rate in the steady-state, frequency, and time domains. In order to compare this to Monte Carlo simulations or experiments the reflectance and the transmittance have to be calculated. For the reflectance in the steady-state domain the fluence rate and its derivative in *z*-direction (flux-term) are used, whereas for all other solutions only the flux term is applied [9, 10, 1].

#### 2.2. Monte Carlo simulations

Monte Carlo simulations are used for the validation of the derived solutions. The basics of our Monte Carlo program were described in the accompanying paper [1] for the calculation of the light propagation in a homogeneous cylinder. For the simulations of the layered cylinder we, additionally, consider the different optical properties of the involved layers using, e.g., Fresnel equations at the boundaries. 10^{7} photons were used in Monte Carlo simulations.

## 3. Results

The derived solutions of the diffusion equation for an N-layered cylinder are first compared to other available solutions of the diffusion equation that are valid for special cases. Then, the derived solutions are verified against results obtained by Monte Carlo simulations.

#### 3.1. Comparison with other solutions of the diffusion theory

In literature so far a solution of the layered cylinder was derived for two and three layers [17, 18]. In Fig. 2 this solution is compared to Eq. 26 for the time resolved reflectance from a finite two-layered cylinder (solid curve). A point source is incident perpendicular onto the center of the cylinder top and the reflectance is calculated at a distance of *ρ* = 17mm.

The disagreement between the time resolved reflectance from the two-layered cylinder calculated with the two different methods is very small, so that we also plot the relative differences, see Fig. 3. These differences are in the range of 10^{-5} and can normally be neglected for applications in the field of biophotonics. If necessary, they can be decreased by putting more effort into the numerical evaluation.

For comparison, also the reflectance from a laterally infinitely extended two-layered turbid medium having the same optical properties as the cylinder is shown in Fig. 2 (dashed curve). As expected, for small times the relative difference is small, whereas for larger times the differences increase due to the fact that the photons remitted at larger times have experienced larger volumes of the turbid media and, thus, have explored the differences of the two considered geometries.

In order to validate the derived solutions for an arbitrary number of layers solutions for different geometries have to be used, because no other solutions for an N-layered cylinder are known so far. Thus, we used the derived solutions of the N-layered turbid cylinder for an infinite large radius (*a* → ∞) and compared it to solutions derived for N-layered turbid media that are laterally infinitely extended [9, 10]. Fig. 4 shows the time resolved transmittance from a 7-layered turbid medium at three different detector distances from the middle of the cylinder bottom (*ρ* = 50, 60, 70mm). The point source is incident at the center of the cylinder top. The optical properties and thicknesses of the seven layers used for the calculations shown in Fig. 4 are listed in Table 1.

Tab. 1: Reduced scattering and absorption coefficients, refractive index, and thickness of the 7-layered turbid medium.

layer |
μ′_{s}[mm^{−1}] |
μ[mm_{a}^{−1}] | n |
l
_{1}[mm] |
---|---|---|---|---|

1 | 1.2 | 0.01 | 1.4 | 3 |

2 | 1.1 | 0.03 | 1.3 | 2 |

3 | 1.3 | 0.02 | 1.5 | 2 |

4 | 1.6 | 0.015 | 1.6 | 2 |

5 | 1.5 | 0.008 | 1.1 | 2 |

6 | 1.4 | 0.035 | 1.7 | 2 |

7 | 1.7 | 0.025 | 1.4 | 3 |

The agreement between the two solutions for all three distances is again excellent. Fig. 5 shows also the relative difference of the two solutions.

For the most time values the relative differences are in the range of 10^{-10} and 10^{-8}. Finally, we compare the derived solutions of the two-layered diffusion equation using the same optical properties in both layers to the solution of a homogeneous cylinder [1]. The point source is incident at the center of the cylinder and time resolved transmittance from the cylinder barrel at a depth of *z* = 7mm is calculated, see Fig. 6. The results are given for two radii of the cylinder (*a* = 18mm (red curve), *a* = 40mm (green curve)). As in the other comparisons the relative difference between the curves calculated with the two models depend on the absolute reflectance data. In general, the larger the reflectance the smaller the relative differences.

The agreement is very good so that we also show the relative difference, see Fig. 7. The red curve shows the differences for *a* = 18mm and the green curve those for *a* = 40mm.

For comparison purpose the transmittance data obtained by a homogeneous parallelepiped model [21] having the same optical properties as the cylinder and a length of 36mm or *a* = 80mm in *x*- and *y*-direction and a length of 10mm in *z*-direction are also shown in Fig. 6 (dashed curves). Similar as for the cylinder the parallelepiped is illuminated at the center and the transmittance is calculated at a depth of *z* = 7mm in the middle of a lateral side, so that the distance between the source and the detector is the same as for the cylinder. For small times the solutions are close, for larger times they again diverge due to the different geometries.

#### 3.2. Comparison with Monte Carlo simulations

In order to test the applicability of the derived solutions for evaluating experimental data we compare them to Monte Carlo simulations. First, the spatially resolved transmittance from a 7-layered cylinder which is illuminated at the center of the cylinder top is investigated. The optical and geometrical properties are those given in Tab. 1, but the refractive index is now *n* = 1.0 for all layers and the surrounding. The transmittance from the center of the cylinder bottom is calculated for three different radii *a* = 4mm (red curve), *a* = 5mm (green curve) and *a* = 10mm (blue curve). Fig. 8 shows the results obtained from the solution of the diffusion equation (solid curves) and those obtained from the Monte Carlo simulations (points).

Good agreement is found for all three radii. This is expected because the distances between the incident source and the detector positions are large so that the diffusion approximation is valid.

Finally, the solution of the diffusion equation is compared to Monte Carlo simulations for a point beam which is non-centric incident onto the cylinder top at *r⃗*
_{0} = (*x*,*y*,*z*) = (4√2mm, 4√2mm, 1/(*μ*
_{a1} + *μ*′_{s1})mm). Fig. 9 shows the spatially resolved reflectance along the *x*-axis (blue curve), along the negative *y*-axis (green curve), and along the bisecting line between the *x*-axis and *y*-axis (red line), see inset. The optical and geometrical properties are those used for the calculations shown in Fig. 8 besides the radius of the cylinder, which is changed to *a* = 12mm.

The spatially resolved reflectance along the bisecting line has a maximum at a distance of 8mm which corresponds to the location of the incident point source, where the diffusion equation fails to describe correctly the reflectance. At large distances, close to 12mm the reflectance from all cylinders rapidly decreases due to the influence of the cylinder border. In general, the agreement between diffusion theory and Monte Carlo simulations is in the expected range.

## 4. Discussion

In this study the solutions of the diffusion equation for an N-layered cylinder were derived in the steady-state, frequency, and time domains applying formulae that we presented in the accompanying paper [1]. Solutions were obtained for a point source incident onto an arbitrary location of the first layer and a circular flat beam incident onto the middle of the cylinder top. In the time domain we derived fast solutions (≈ 10ms) for the special case of equal reduced scattering coefficient and refractive index in each layer. For the general case in the time domain the calculation time is about 1s when using a state-of-the-art computer and the Pascal (Delphi) programming language, whereas in the steady-state domain the calculation times are orders of magnitude faster.

A solution up to three layers was reported in the literature using a different approach [18]. We compared this solution to the one derived in this paper for two-layers (and also three-layers, data not shown) and found excellent agreement. We also compared the derived solutions for the special cases of infinitely large radii of a 7-layered turbid medium and of homogeneous cylinders and observed very good accordance.

For the validation against transport theory we compared the derived solutions with Monte Carlo simulations. Calculations in transmission and remission were performed showing the expected small relative differences in the range of a few percent.

In general, although the derivation of the above presented solutions of the diffusion equation might be somewhat cumbersome, the obtained final results can be easily implemented. In addition, we provide executable files of the solutions in the steady-state and time domains on the internet [22].

In the future, the solutions will be used to examine the inverse problem, the determination of the optical properties (and layer thickness) of (some of) the layers from measurements in the different domains.

## Acknowledgement

We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076). We like to thank Fabrizio Martelli for providing us his code to solve the two-layered diffusion equation used for the comparison in Fig. 2.

## References and links

**1. **A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

**2. **A. Ishimaru, “Wave Propagation and Scattering in Random Media,” Academic Press, New York (1978).

**3. **I. Dayan, S. Havlin, and G.H. Weiss, “Photon Migration in a Two-Layer Turbid Medium. A Diffusion Analysis,” J. Mod. Opt. **39**, 1567–1582 (1992). [CrossRef]

**4. **A. Kienle, M.S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. **37**, 779–791 (1998). [CrossRef]

**5. **A. Kienle, T. Glanzmann, G. Wagnières, and H. van den Bergh, “Investigation of Two-Layered Turbid Media with Time-Resolved Reflectance,” Appl. Opt. **37**, 6852–6862 (1998). [CrossRef]

**6. **A. Kienle and T. Glanzmann, “In Vivo Determination of the Optical Properties of Muscle Using a Layered-Model,” Phys. Med. Biol. **44**, 2689–2702 (1999). [CrossRef] [PubMed]

**7. **J.M. Tualle, H.M. Nghiem, D. Ettori, R. Sablong, E. Tinet, and S. Avrillier, “Asymptotic Behavior and Inverse Problem in Layered Scattering Media,” J. Opt. Soc. Am. A **21**, 24–34 (2004). [CrossRef]

**8. **X.C. Wang and S.M. Wang, “Light Transport Modell in a N-Layered Mismatched Tissue,” Waves Rand. Compl. Media **16**, 121–135 (2006). [CrossRef]

**9. **A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt., accepted. [PubMed]

**10. **A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt., accepted. [PubMed]

**11. **G. Alexandrakis, T.J. Farrell, and M.S. Patterson, “Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium,” Appl. Opt. **37**, 7401–7409 (1998). [CrossRef]

**12. **S-H. Tseng, C. Hayakawa, B.J. Tromberg, J. Spanier, and A.J. Durkin, “Quantitative Spectroscopy of Superficial Turbid Media,” Opt. Lett. **23**, 3165–3167 (2005). [CrossRef]

**13. **G. Alexandrakis, T.J. Farrell, and M.S. Patterson, “Monte Carlo Diffusion Hybrid Model for Photon Migration in a Two-Layer Turbid Medium in the Frequency Domain,” Appl. Opt. **39**, 2235–2244 (2000). [CrossRef]

**14. **M. Das, C. Xu, and Q. Zhu, “Analytical Solution for Light Propagation in a Two-Layer Tissue Structure with a Tilted Interface for Breast Imaging,” Appl. Opt. **45**, 5027–5036 (2006). [CrossRef] [PubMed]

**15. **A.H. Barnett, “A Fast Numerical Method for Time-Resolved Photon Diffusion in General Stratified Turbid Media,” J. Comp. Phys. **201**, 771–797 (2004). [CrossRef]

**16. **C. Donner and H.W. Jensen, “Rapid Simulations of Steady-State Spatially Resolved Reflectance and Transmittance Profiles of Multilayered Turbid Materials,” J. Opt. Soc. Am. A **23**, 1382–1390 (2006). [CrossRef]

**17. **F. Martelli, A. Sassaroli, S. Del Bianco, Y. Yamada Y, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for Layered Diffusive Media by the Eigenfunction Method,” Phys. Rev. E **67**, 056623 (2003). [CrossRef]

**18. **F. Martelli, A. Sassaroli, S. Del Bianco, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for a Three-Layer Medium: Application to Study Photon Migration through a Simplified Adult Head Model,” Phys. Med. Biol. **52**, 2827–2843 (2007). [CrossRef] [PubMed]

**19. **F. Martelli, S. Del Bianco, and G. Zaccanti, “Perturbation Model for Light Propagation through Diffusive Layered Media,” Phys. Med. Biol. **50**, 2159–2166 (2005). [CrossRef] [PubMed]

**20. **F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, “Light Propagation through Biological Tissue and other Diffusive Media,” SPIE Press, Bellingham (2010).

**21. **A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A **22**, 1883–1888 (2005). [CrossRef]