## Abstract

The efficiency of using intensity modulated light for the estimation of scattering properties of a turbid medium and for ballistic photon discrimination is theoretically quantified in this article. Using the diffusion model for modulated photon transport and considering a noisy quadrature demodulation scheme, the minimum-variance bounds on estimation of parameters of interest are analytically derived and analyzed. The existence of a variance-minimizing optimal modulation frequency is shown and its evolution with the properties of the intervening medium is derived and studied. Furthermore, a metric is defined to quantify the efficiency of ballistic photon filtering which may be sought when imaging through turbid media. The analytical derivation of this metric shows that the minimum modulation frequency required to attain significant ballistic discrimination depends only on the reduced scattering coefficient of the medium in a linear fashion for a highly scattering medium.

© 2016 Optical Society of America

## 1. Introduction

Imaging through and within turbid media is an area of interest that has tremendous application in medical diagnostics [1], underwater vision [2], imaging through colloids [3] and transportation and navigational aids [4, 5]. Light traveling through a complex medium with randomly distributed positions and refractive indices undergoes absorption and random scattering and loses the spatial and temporal information of its source. The photons that undergo such multiple scattering are labeled as diffused photons. A small fraction of the total photons called ballistic and quasi-ballistic photons, undergo no or very few forward scattering events before they reach the detector, and they retain the information of its source (direction, polarization state, modulation,…). It is of wide interest to discriminate the ballistic/quasi-ballistic photons from the diffused photons for resolution enhanced imaging through turbid media. However, the diffuse light that strongly depends on the properties of the scattering medium can be used to deduce various parameters related to the medium itself. Thus, imaging in turbid media can be classified into two broad categories: parameter estimation using only diffused light to obtain image of heterogeneities in the turbid media and filtering of ballistic photons from diffused light for high resolution imaging through turbid media.

#### Parameter estimation

In parameter estimation, the physical properties of the intervening media are of interest and controlled light sources may be used to estimate the scattering and absorption parameters [6, 7]. Diffuse optical imaging which has been widely studied and applied in medical imaging [8] is used to estimate the scattering and absorption coefficient in the intervening media for detection of malignant tissues in breast [9] or for brain imaging [10]. Various methodologies have been developed that utilize the diffusion equation to model light transport through tissues to estimate its scattering parameters. Steady-state solution is generally used in case where continuous wave light sources are used [11] or when patterned light source is used as in case of spatial frequency domain imaging [12]. Diffusion approximation is also commonly used to model the transport of pulsed light [13, 14] or modulated light [15–18] through the scattering media. The precision of estimation of the parameters is crucial in this case. In practice, intensity modulated light with diffusion theory are widely used in frequency-domain photon migration measurements, where an intensity modulated light with modulation frequencies of a few hundred MHz is transmitted through a diffusing medium. According to the time dependent solution of the diffusion theory, at these frequencies, the wave number of the density wave traveling through the medium is dependent on the scattering and absorption properties of the medium. As a result, the detected modulated light has reduced modulation index and an additional phase, both of which depend on the scattering properties of the medium and its thickness. As the diffusion theory provides an analytical model for the change of modulation index and phase of the modulated light, the parameters of the intervening medium could be estimated by using a single or multiple modulation frequencies.

#### Ballistic discrimination

On the other hand, ballistic filtering is used when the spatial resolution of objects embedded in scattering media is of interest [19, 20]. For instance, for industrial applications where it is required to image objects embedded in colloidal system, or in navigation where vision through fog can be efficiently achieved with ballistic filtering. The problem of imaging through such media has been addressed using various techniques that essentially rely on discriminating ballistic/quasi-ballistic photons from multiply scattered (diffuse) photons traversing through a turbid medium. For example, time gated imaging [21,22], polarization gated imaging [23], intensity modulation imaging [24] etc., where information carried by the ballistic/quasi-ballistic photons is filtered from the diffuse light. It is challenging to efficiently detect the ballistic photons that have low signal to noise ratio as compared to the diffused photons, especially in highly scattering media. Ballistic photon filtering could also be achieved with an intensity modulation imaging scheme by discriminating them from the diffused light that arrives with reduced modulation index and additional phase. One of the purposes of this article is to theoretically study the feasibility of such an approach, which is currently the topic of active research [25].

In this article, we analyze the widely used diffusion approximation for transport of modulated light, from an information theoretical point of view for efficient parameter estimation and for ballistic photons discrimination. Assuming a general scheme of optical quadrature demodulation detection and a corresponding noise model, this approach provides a rigorous insight into the maximal precision of estimation of parameters using the diffusion equation, thereby providing a robust theoretical argument for choosing the frequency of modulation suited to the experiment at hand. Similarly, for ballistic filtering, we will present an information theoretic performance metric and study the optimum ballistic discrimination efficiency that can be achieved under this imaging scheme.

In section 2, we provide a brief description of the use of diffusion equation, the transport of modulated light in diffusing media and the imaging scheme along with the parameters that will be used in the theoretical analysis. In section 3, a generalized demodulation scheme is defined and a noise model is presented. Then, in section 4, the noise model and the diffusion theory are used to calculate the Cramer-Rao bound for parameter estimation. Finally, in section 5, the ballistic filtering efficiency is analyzed and discussed with respect to real field situations.

## 2. Imaging scheme and diffusion model

The diffusion theory for photon transport provides a simple, fast and analytical method for modeling the light propagation through various turbid media. The properties of intensity modulated light through a diffusing medium have been well studied and reviewed [15, 26–28]. Here, we present a brief introduction underlining concepts that are relevant in the context of this article. We will follow the modeling of the intensity modulated light in diffusing media as derived in several references [16, 28, 29] in the context of diffuse optical imaging. We will use the formulations previously derived by the cited authors to parameterize our detection scheme and then use information theoretical tools to derive simple rules for choosing frequency operating points to efficiently use intensity modulation light depending on the properties of the medium.

#### 2.1. Diffusion model

The diffusion model arises when the photons are allowed to perform a random walk, diffusing from high photon density regions to low photon density regions. The theory has proved efficient when modeling light in a predominantly scattering medium where the source and the detector are far from boundaries of the medium and when detection is carried out sufficiently away from a point source, in a supposed infinite geometry. The efficiency of the diffusion theory has been studied alongside Monte-Carlo simulation and shown to have well acceptable accuracy in the domain of its validity [30]. Let us now identify some parameters that are important for describing light transport in a scattering medium. First, the absorption coefficient, denoted by *µ* (in units of *m*^{−1}), corresponds to the inverse of the mean distance traveled by a photon before it is absorbed. Similarly, the scattering coefficient *µ _{s}* (

*m*

^{−1}) is the inverse mean distance before a photon is scattered in the medium. We shall also consider the

*reduced*scattering coefficient, denoted by

*σ*in this work, and classically defined as

*σ*=

*µ*(1−

_{s}*g*), where

*g*is the (dimensionless) anisotropy factor, defined as the mean cosine of the scattering phase function [31,32]. Such a reduced scattering coefficient allows anisotropic scattering media to be easily described (notably through diffusion equation), and may be interpreted as the inverse of the scattering length for a fictitious isotropically scattering (

*g*= 0) medium having the same macroscopic light-transport parameters as the real anisotropically scattering (

*g*≠ 0) medium. Within this framework, an important length scale is the transport mean free path (

*TMFP*) denoted by

*l*, which is the mean distance traveled by photons before they lose their initial directional information and is equal to [

_{*}*µ*+

*σ*]

^{−1}. The diffusion theory also includes two other constants: the diffusion length (

*D*) defined as

*D*=

*l*

^{*}/3 and the optical penetration depth (

*δ*), which is the inverse of the effective attenuation coefficient $\left(\sqrt{\mu /D}\right)$ in diffusing medium.

In addition, for the analysis in the following sections we will use various dimensionless constants that are also given in the right side of Table 1. For instance, the dimensionless parameter *R _{δ}* =

*r/δ*corresponds to the effective optical attenuation of diffused light. The parameter

*q*is related to the angular frequency of modulation and the non-trivial form of the reparameterization will be justified in following sections. For ease of reading, the corresponding definitions of all these parameters are tabulated in Table 1.

#### 2.2. Imaging scheme

The imaging scheme considered in this article includes a directional point source of light with a forward cone solid angle Ω and that subtends a solid angle of *d*Ω at the detector which is placed at a distance *r* from the source. For the sake of simplicity, we shall consider a source of limited spectral range *λ*_{0} *±* Δ*λ*, so that the above diffusion parameters can be considered as constant over the considered spectral range. The schematic in Fig. 1 illustrates the scenario. In the presence of an intervening scattering medium, the net intensity of ballistic photons reaching the detector of collection area *d*Ω*r*^{2} is proportional to the total power *P*_{0} emitted by the source. It is given by the Bouguer-Beer-Lambert’s law as
${I}_{B}=\xi {P}_{0}{e}^{-(\mu +{\mu}_{s})r}\frac{d\mathrm{\Omega}{r}^{2}}{\mathrm{\Omega}{r}^{2}}=\xi {P}_{0}{e}^{-(\mu +{\mu}_{s})r}\frac{d\mathrm{\Omega}}{\mathrm{\Omega}}$ [31], with an extinction coefficient *µ* + *µ _{s}*, (equal to the inverse MFP, see Table 1) allowing us to define

*R*=

_{b}*r/MFP*as the effective optical attenuation of ballistic light. The scaling factor

*ξ*represents the overall detector efficiency on the considered spectral range. Similarly, according to the steady state solution of the diffusion theory, the diffuse photon intensity reaching the detector is also proportional to the total power emitted by source, such that ${I}_{D}=\xi {P}_{0}{e}^{-r/\delta}\frac{d\mathrm{\Omega}{r}^{2}}{4\pi Dr}=\xi {P}_{0}{e}^{-r/\delta}\frac{d\mathrm{\Omega}}{4\pi}\frac{r}{D}$ [16, 29]. The expressions for the intensities of ballistic and diffused light show that there are clearly two important length scales to be considered, namely, the TMFP (

*l*

^{*}) and the optical penetration depth (

*δ*). Considering only the above two classes of photons, it is possible to obtain an order of magnitude value of the ratio (

*α*) of ballistic photons to diffuse photons reaching the detector as

*R*

^{*},

*R*,

_{δ}*R*and Ω′ = 4

_{b}*π/*Ω.

#### 2.3. Modulated light in diffusing media

The propagation of sinusoidally modulated light through a scattering medium has been modeled using diffusion theory and it has been shown that the transport of modulated light behaves as photon density waves whose properties are dependent on the properties of the medium [15,16, 27, 28]. In this article, we consider an intensity modulated source of light having modulation angular frequency *ω* and amplitude modulation index *M* (also classically termed *modulation depth*), which describes the amplitude of the periodic fluctuation of light intensity around its unmodulated level. In that case, the instantaneous intensity reads *i*(*t*) = *I*_{0}(1 + *M* cos[*ωt*]). The ballistic light that follows Bouguer-Beer-Lambert’s law, is only attenuated and reaches the detector with instantaneous intensity *i _{b}*(

*t*) =

*I*(1+

_{B}*m*cos[

_{B}*ωt*]) without any change in received modulation index,

*m*=

_{B}*M*.

However, the time dependent solution of the photon diffusion theory shows that the modulated light traversing through a scattering medium is received at the detector with reduced modulation index and additional phase [16]. Then, the instantaneous diffuse light intensity received at the detector is *i _{d}*(

*t*) =

*I*(1 +

_{D}*m*cos[

_{D}*ωt*+ Δ

*φ*]). Without derivation, we present the expression of the reduced modulation index

*m*and the phase Δ

_{D}*φ*which is identically reported in [16,29] for an infinite medium:

*ω*ϵ [0,∞). Although the physical interpretation of this parameter is not straightforward, we will see later that

*qR*can be identified as a dimensionless, frequency-dependent, effective attenuation of the diffused light. In the remainder of this article, we define

_{δ}*β*=

*m*as the ratio of the modulation indices of ballistic light to diffuse light. The expressions of the intensity, phase and modulation index of ballistic and diffuse light components are recalled in Table 2.

_{B}/m_{D}It is quite straightforward to see that using the above equations, the parameters *R _{δ}* and

*q*can be estimated when the modulation index and phase of the diffuse light are accurately detected. The model can indeed be inverted as shown in Eq. (3):

Thus, the above formulation provides a simple analytical method for estimation of the scattering and absorption parameters of the scattering medium using only diffuse light and a modulated light source. In practice, especially in diffuse optical imaging, the modulation frequency is scanned to obtain corresponding values of modulation index and phase. Then, a non-linear fit of the theoretical prediction with the data provides the estimates for the scattering and absorption properties of the medium [16]. The effect of scattering media on modulation index and phase can also be exploited to attain discrimination of ballistic photons that retain the modulation index and phase. These application scenarios can be analyzed from an information theoretical point of view for a well-defined detection technique. The detection of the modulation index and the phase can be performed in various ways. Generally in a demodulation scheme, the amplitude, phase and mean intensity are recorded and then the modulation index can be easily estimated. Quadrature demodulation is one of the simplest and the most widely used scheme for demodulation. It avoids, in particular, phase tracking of the incoming signal which brings additional noise contributions. In the following section we look at a quadrature detection scheme and derive the noise model for the detection.

## 3. Quadrature detection scheme and noise model

Demodulation of a sinusoidally modulated signal can be performed by product detection, where the signal is continuously multiplied with a sine and cosine of the same frequency to obtain the quadrature amplitudes. This quadrature detection can be interpreted as follows. Let us consider a detection system where a sinusoidal signal of angular frequency *ω*, modulation index *M* and mean intensity *I*_{0} is sampled *N* times in a period producing photon count or grey level data *d _{i}* (

*i*= 1 to

*N*). Then the quadrature components can be estimated by applying sinusoidal weighting of the data points and computing the two quadratures $U={\displaystyle {\sum}_{i=1}^{N}{d}_{i}}\mathrm{cos}\left[\omega {t}_{i}+\varphi \right]=A\mathrm{cos}\left[\varphi \right]$ and $V={\displaystyle {\sum}_{i=1}^{N}{d}_{i}}\mathrm{sin}\left[\omega {t}_{i}+\varphi \right]=A\mathrm{sin}\left[\varphi \right]$, for

*N ≥*2, where

*ϕ*denotes the relative phase difference between the signal and the demodulating waveforms. From the obtained variables

*U*and

*V*, estimates of the amplitude

*A*=

*MI*

_{0}/2 and relative phase

*ϕ*of the demodulated signal are respectively given by $\widehat{A}=\sqrt{{U}^{2}+{V}^{2}}$ and $\widehat{\varphi}={\mathrm{tan}}^{-1}\left[V/U\right]$.

The intensity statistics of the random variables *U* and *V* can be derived by knowing that the photon count within an infinitesimal sampling window at time *t _{i}* may follow Poisson distribution with mean and variance equal to the intensity received at the sampling window, which is

*d*in the current notation. Then, a weighted addition of the Poisson distributed data is also Poissonian, indicating that

_{i}*U*and

*V*will also be distributed as Poisson random variables with variance

*I*

_{0}/2 (see Appendix A.1). For high intensities, the Poisson distribution can be well approximated by a Gaussian distribution $\left(\mathcal{N}\right)$ with variance equal to the mean. Thus, we approximate that

*U*and

*V*are distributed as $\mathcal{N}\left(A\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\left[\varphi \right],{\mathrm{\Lambda}}^{2}\right)$ and $\mathcal{N}\left(A\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left[\varphi \right],{\mathrm{\Lambda}}^{2}\right)$, respectively, with Λ

^{2}=

*I*

_{0}/2 and $\mathcal{N}\left(\overline{x},var\left(x\right)\right)$ denoting the normal distribution with mean $\overline{x}$ and variance

*var*(

*x*). Knowing the distribution of the quadrature components, it is then possible to derive the joint probability density function (PDF) of the random variables, $Z=\sqrt{{U}^{2}+{V}^{2}}$ and Ψ = tan

^{−}^{1}(

*V/U*), respectively associated with the modulation amplitude

*A*and relative phase

*ϕ*, by applying the appropriate change of variables to the Gaussian distribution, as given in the Appendix A.2. The expected values of the random variables

*Z*, Ψ and their variance form the parameter vector

*θ*′ = [

*A,ϕ*,Λ

^{2}], where 2

*A*=

*MI*

_{0}is the peak-to-peak amplitude as shown in the schematic of Fig. 2. Then, the joint PDF of the amplitude and phase is

Quadrature demodulation cameras are in widespread use, especially in 3D imaging. For instance, novel demodulation cameras like time-of-flight (TOF) cameras collect four samples per period to obtain the quadratures *U* = *d*_{1} *−d*_{3} and *V* = *d*_{2} *−d*_{4} at frequencies up to 20 MHz [33]. Information theoretical studies have been made to accurately determine the phase which contains the depth information of a scene using the above PDF [34]. In the following, we will use Fisher information (FI) and Cramer-Rao bound (CRB) to analyze how well the diffusion parameters of the scattering medium can be estimated using modulated light and the quadrature detection technique.

#### 3.1. Fisher information

Fisher information (FI) is a useful and efficient method of quantifying the precision with which the unknown parameters in a data model can be estimated. It can be further used to analyze the behavior of the covariance of estimators and its dependency on other parameters. FI is ultimately related to the minimum covariance bound on the unbiased estimation of the unknown parameters. Here, we will use the FI to quantitatively derive simple rules that must be taken into account when using intensity modulated light in diffusion approximation for estimation of scattering properties of a medium and for ballistic filtering applications.

The expected Fisher information matrix (FIM) for a parameter vector (*η*) is defined as
$\mathcal{F}{\left(\eta \right)}_{ij}=-\u3008\frac{{\partial}^{2}\mathrm{ln}[P(\mathbf{x}|\eta )]}{\partial {\eta}_{i}\partial {\eta}_{j}}\u3009$, where 〈.〉 denotes the expectation value. We present below the FIM for the detection procedure described above with respect to the parameter vector *θ*′ = [*A, ϕ*,Λ^{2}]. It is straightforward to show that the FI reads in that case

The above FIM is diagonal, which shows that the three parameters in *θ*′ can be estimated independently. Moreover, the precision in estimation of the phase increases with the amplitude of the signal. The detection technique described in the preceding section is not limited to detector arrays, but is generally adopted in lock-in detection. In this detection scheme, the recorded amplitude and phase of the diffuse light through a scattering media can be modeled by the diffusion theory as presented in preceding sections. Using this model, the noise model of the detection can be reparametrized as shown below and consequently, estimation precision of the parameters can be calculated.

#### 3.2. Reparametrization of noise model using diffusion theory

Let us consider the effect of diffused light only, and obtain the FIM with respect to the new set of parameters *θ* = [*M,R _{δ},R*

^{*}] to deduce some insight into the precision of estimation of each respective parameter and its dependency on other parameters. The FI with respect to the new parameters can be obtained by calculating the Jacobian matrix, ${I}_{\mathcal{D}}$, of the transformation

*θ*′

*→ θ*in the presence of diffuse light only, as denoted by the subscript $\mathcal{D}$. The

*i*component of this matrix is ${\left[{J}_{\mathcal{D}}\right]}_{i,j}=\frac{\partial {{\theta}^{\prime}}_{i}}{\partial {\theta}_{j}}$. Given the modulation index, phase and the intensity expected for diffuse light from the diffusion model, the amplitude, phase and variance can be written as ${A}_{D}={m}_{D}{I}_{D}/2,\mathrm{\Delta}\varphi ={R}_{\delta}\sqrt{2\left({q}^{2}-1\right)}$ and ${\mathrm{\Lambda}}_{D}^{2}={I}_{D}/2=3{S}_{0}{R}_{*}{e}^{-{R}_{\delta}}/2$, where we have set

^{th}, j^{th}*S*

_{0}=

*ξd*Ω

*P*

_{0}/4

*π*. The Jacobian is then calculated to obtain the FIM with respect to parameter vector

*θ*using the transformation ${\mathcal{F}}_{\mathcal{D}}(\theta )={J}_{\mathcal{D}}^{T}\mathcal{F}({\theta}^{\prime}){J}_{\mathcal{D}}$. Both, the Jacobian and the FIM ${\mathcal{F}}_{\mathcal{D}}(\theta )$ are shown in Appendices B and C, for reference. It is worth mentioning here that in applications where the SNR of the amplitude and phase are important, the noise model presented here can be used to obtain the SNR of detected amplitude and phase as a function of frequency by using the above relations. However, we will focus on deriving and analyzing the minimum variance bounds on estimation of scattering parameters as presented in the following section.

## 4. Maximal precision in estimation of scattering parameters in diffuse optical imaging

#### 4.1. Lower bound on estimation variance

The above re-parameterization allows us to study the variance bound in estimation of *θ* with respect to the frequency of modulation represented by the variable *q*. According to the Cramer-Rao theorem, for any unbiased estimator
$\widehat{\theta}$ of the parameter vector *θ*, and for any row vector *w* having the same dimension as *θ*, one has
$w\u3008\widehat{\theta}{\widehat{\theta}}^{T}\u3009{w}^{T}\ge w\mathcal{F}{(\theta )}^{-1}{w}^{T}$. Thus, the Cramer-Rao bound (CRB) provides a minimum covariance bound that can be reached by an efficient estimation technique. Generically, it is not guaranteed that such an efficient estimator always exists, however, Maximum Likelihood (ML) estimators has been shown to be asymptotically unbiased and efficient for a large collection of data, with optimality strictly proved in the case of exponential family of distributions [35]. The exact forms of ML estimators in this case are solutions to transcendental equations and remain out of the scope of this article. Instead, let us analyze the lower covariance bound to reveal some simple conclusions and rules for the estimation problem. We provide the bound as
$CR{B}_{\mathcal{D}}\left(\theta \right)={\mathcal{F}}_{\mathcal{D}}{\left(\theta \right)}^{-1}$ for the invertible matrix
${\mathcal{F}}_{\mathcal{D}}\left(\theta \right)$ which is shown in the Appendix C:

It is clearly seen that the variances are functions of the frequency of modulation (represented by *q*). Let us first look at the variance in estimation for the parameters R* _{δ}* and R

_{*}. The variance of

*R*decreases with

_{δ}*R*

^{*}. The dependence of both the parameters on frequency is the same and has a functional form ${e}^{\left(-1+2q\right){R}_{\delta}}/\left(-1\phantom{\rule{0.2em}{0ex}}+{q}^{2}\right)$. A simple calculus shows that this function reaches a minimum at ${q}_{opt}=(1+\sqrt{1+4{R}_{\delta}^{2}})/2{R}_{\delta}$. This indicates that there exists an optimal angular frequency

*ω*around which the variance of the estimation is minimum. This optimal angular frequency depends only on

_{opt}*R*and is given by the following expression

_{δ}*R*limit, the optimal frequency behaves like ${\omega}_{opt}/\mu c=2/\sqrt{{R}_{\delta}^{2}}$, whereas it evolves as ${\omega}_{opt}/\mu c=2/\sqrt{{R}_{\delta}}$ when

_{δ}*R*∞. In the next subsection we analyze the properties and evolution of the optimal frequency

_{δ}→*ω*obtained for the estimation of scattering parameters.

_{opt}#### 4.2. Optimal operating frequency

The expression for the variance-minimizing optimal angular frequency *ω _{opt}* has a non-trivial form that depends on the normalized optical attenuation

*R*=

_{δ}*r/δ*, which basically represents the inverse of an overall visibility factor for the diffuse photons. The existence of an optimal operating point has interesting consequences, especially in the context of diffused optical imaging, where the SNR of the images is an important limiting factor. Generally, in such applications, the operating frequency, denoted here by

*ω*, is chosen such that

_{a}*ω*= 1 [29]. This is indeed a reasonable choice of operating frequency, justified by the fact that at smaller frequencies the phase change due to diffusion is too small to be detected, while at much higher frequencies, the change in phase becomes comparatively insensitive to further increase in frequency, thereby not bringing any further improvement in the estimation of the medium parameters [29].

_{a}/µcTo understand the loss incurred by an imperfect choice of operating frequency, we compare the performance of estimation provided by the two frequencies, *ω _{opt}* and

*ω*. We plot the ratio of the variance in the estimation of

_{a}*R*obtained when using

_{δ}*ω*against

_{opt}*ω*in Fig. 3(b). The figure provides a quantitative study of the loss in optimal precision in using

_{a}*ω*as opposed to using

_{a}*ω*. The precision of estimation are equal only when

_{opt}*R*= 5.3, at which point

_{δ}*ω*=

_{opt}*ω*, as indicated by the dashed lines in the figures. For example, taking typical values of scattering parameters valid in tissues, such that

_{a}*µ*= 0.1 cm

^{−1}and

*σ*= 10 cm

^{−}^{1}, the two operating points are the same only for a detector placed as distance of 3 cm. The standard operating frequency

*ω*is independent of the propagation distance

_{a}*r*and of

*σ*, depending only on the absorption coefficient

*µ*. As illustrated in Fig. 3(a), for a detection carried out at a distance higher than 3 cm, the optimal frequency is in fact smaller, which may reduce the cost and complexity of the electronics required at high frequencies while providing better results. On the other hand, if the detection is carried out closer than 3 cm, operating frequencies higher than the standard operational frequency are suggested.

Similarly, since *R _{δ}* is a function of

*σ*,

*µ*and

*r*, the dependency of the optimal modulation frequency on these three parameters can be analyzed from the contour plot shown in Fig. 3(c). Firstly, for fixed values of

*σ*and

*µ*, the optimal operating frequency decreases with increase in

*r*. Secondly, for fixed values of

*r*and

*µ*, the optimal frequency is seen to decrease with the increase in

*σ*. The above two dependencies can be interpreted by noticing that longer distance of travel or higher scattering coefficient would allow for a larger amount of multiple scattering events leading to greater modification of the phase and the modulation index of the diffused light, for a given coefficient

*µ*.

Finally, the optimal frequency is seen to increase with the increase in *µ*, for fixed values of *r* and *σ*. An increase in absorption would lead to smaller amplitude of modulation at the detector. To compensate for this loss that occurs at a rate of *µc*, the photon density arrival rate should be increased leading to a increase in optimal frequency of operation. It can also be noticed that the change of optimal operating frequency is more sensitive to a change in absorption coefficient than to a change in *σ* or *r*. The above conclusions can also be easily obtained by noticing that equation Eq. (7) can be approximated to
${\omega}_{opt}=2\mu c/\sqrt{{R}_{\delta}}=2\sqrt{\mu}c{\left[r\sqrt{3\left(1+\sigma /\mu \right)}\right]}^{-1/2}$, for sufficiently large value of *R _{δ}*, as noted above. It is interesting to notice at this level that the above analytical deductions about the existence of an optimal operating frequency and its evolution with scattering parameters are qualitatively in very good agreement with previous numerical simulations and experimental studies presented in [36–38].

More generally, the above results show that our analysis, based on a diffusion model for the propagation of photons coupled with an information theoretical approach, makes it possible to account for several competing phenomena involved in light diffusion. For instance, the noise variance Λ^{2} observed with the quadrature detection scheme described in Section 3 decreases with *R _{δ}* but is independent of

*ω*, whereas the modulation index and the phase depend on

*R*but with rates of change that are functions of the angular frequency). This analysis is able to provide the functional dependence of the optimal operating frequency with respect to properties of the medium under consideration and practical indications concerning the optimal experimental parameters for the estimation task considered. The conclusion obtained appears to be more specific than the usual rule of thumb used in similar experiments, as it reveals the optimal operating frequency and provides insight into the various competing phenomena that produce it.

_{δ}It is possible to extend such an information theoretical analysis towards another problem, namely, ballistic photon filtering in diffusing media, where the goal is to efficiently discriminate the ballistic photons from the diffused photons. In the next section we will introduce a performance metric for ballistic discrimination and present the optimal operating point based on the information theoretical approach coherent with the previous discussions.

## 5. Ballistic filtering

The importance of ballistic discrimination has been discussed in the introduction section, for instance when imaging objects embedded in nebulous media like fog. As noted above, when high-frequency modulated light propagates in a turbid medium, the ballistic photons and diffused photons have different transport properties in terms of modulation index and phase, which can be exploited by a demodulating detector to attain discrimination of ballistic photons.

By contrast, such a physical ballistic filtering effect cannot be obtained with a standard intensity camera, or with low-frequency modulation/lock-in detection. In both cases, detection of an extremely small ballistic contribution over a spatially uniform diffuse illumination would be possible only at the expense of a dramatic increase in the detector dynamics or acquisition time.

In this section, we investigate the conditions required to obtain significant ballistic filtering with modulated light in a turbid medium. More particularly, we derive the minimum modulation frequency required to attain ballistic filtering irrespective of photon budget and exposure time.

#### 5.1. Gain definition for ballistic filtering efficiency

We will again resort to FI to define a ballistic filtering efficiency as the gain in information provided by the ballistic light for the estimation of the modulation index *M* of the light source. Let us consider a quadrature demodulation camera that receives ballistic light over the diffused light on a set of pixels denoted
$\mathcal{B}\oplus \mathcal{D}$ and another set of pixels that receive only diffused light at a region denoted by
$\mathcal{D}$. In most cases, the contrast between these two regions is marginal because fewer ballistic photons reach the detector. To quantify this contrast in a demodulation scheme, we define the ballistic discrimination efficiency (or gain) as the ratio of FIs in the estimation of the actual modulation index *M* of the source when using data from region
$\mathcal{B}\oplus \mathcal{D}$ as opposed to
$\mathcal{D}$. The gain in information, denoted by
${\mathcal{G}}_{bf}$ is thus defined as

*θ*′

*→ θ*in the presence of ballistic and diffuse light. Here, we do not need to compute either the entire Jacobian ${J}_{\mathcal{B}\oplus \mathcal{D}}$ or the entire FIM ${\mathcal{F}}_{\mathcal{B}\oplus \mathcal{D}}(\theta )$, as our problem is restricted to finding the first term in the FIM, which provides a limiting but reasonable condition for achieving ballistic filtering when other parameters like

*R*and

_{δ}*R*are already known or assumed to be known. Consequently, the gain can be calculated as

_{*}*τ*= 2(

*qR*), $\beta ={e}^{{R}_{\delta}(q-1)}$ and

_{δ}−R_{b}*α*=

*I*is the ratio of received ballistic light over diffused light.

_{B}/I_{D}#### 5.2. Condition for ballistic filtering gain

The contour plot of the
$\mathrm{ln}[{\mathcal{G}}_{bf}]$ as a function of the normalized angular modulation frequency (*ω/µc*) and of *σ/µ* is shown in Fig. 4 for two values of the anisotropy coefficient *g*, taking Ω′ = 1 (isotropic emitter) and *rµ* = 10. In Figs. 4(a) and 4(b), the domain of validity of the diffusion approximation is delimited by the region below the yellow dashed line (see Section 5.4). It can be seen that a gain much greater than unity can be obtained for angular frequencies that lie above a roughly linear contour. Thus, the minimum angular frequency required for ballistic filtering depends roughly linearly on the scattering properties (through *µ _{s}*,

*σ*and/or

*g*) of the medium. The analytical expression of the contour for unity gain is not easily computed because of the oscillating cosine term in Eq. (10a), as can be seen in the inset of Fig. 4. Taking a closer look at the expression of ballistic gain, it can be noticed that a significant gain can be obtained only when

*τ*> 0, which allows exponential increase in gain and at the same time makes it possible to neglect the cosine term inside the bracketed expression. Moreover, it can be easily checked that it is a sufficient condition to ensure that ${\mathcal{G}}_{bf}$ of Eq. (10a) is greater than unity. The above simple condition can be rewritten as

*qR*, which interestingly suggests that

_{δ}> R_{b}*qR*effectively behaves as a frequency-dependent attenuation for diffused light, which should be greater than the effective attenuation experienced by ballistic light (

_{δ}*R*) in order to achieve efficient ballistic filtering. The above inequality leads to a condition in terms of minimum frequency of operation which can be written as

_{b}*r*. When expressed in terms of angular frequency

*ω*, this gain condition can be plotted as the red dashed line on the contour maps in Fig. 4 and seems to provide a well-defined condition for attaining ballistic gain. Further simplification holds under the validity conditions of the diffusion theory, where

*σ*> 10

*µ*. In this case, the condition for achieving ballistic filtering reduces to and is insensitive to the value of the absorption coefficient

*µ*. This clearly indicates that filtering of ballistic photons using an imaging demodulation technique will be more difficult, and hence require higher modulation frequency, in forward anisotropically scattering media with

*g*> 0. Indeed, this can be physically understood, as diffused photons will be less likely to deviate from the ballistic path in a significantly forward scattering medium.

As a result, the above equation provides a simple rule of thumb for achieving ballistic filtering in the context of intensity modulation and quadrature detection. It must be noted that the effect of Ω′ and *rµ* (that are set as constants in Fig. 4) is only on the value of the gain and they do not affect the above condition for efficient ballistic discrimination. Identification of this minimum modulation frequency for ballistic filtering is also important from a practical point of view as it is easier to design electronic or electro-optic devices that work at low modulation and demodulation frequencies. According to the above results, ballistic filtering in biomedical applications would, however, require very high (if not unrealistic) modulation frequencies, as for typical values of *σ* = 10 cm^{−}^{1}, refraction index *n* = 1.33, and a propagation distance *r* = 5 cm, one obtains *f _{min}*=23.9 GHz (respectively

*f*=597 GHz) when

_{min}*g*= 0 (respectively

*g*= 0.8)!

The expression derived above can also serve conversely to provide the range of visibilities that can be handled by a ballistic filtering device working at any fixed modulation-demodulation frequency. For instance, for transport safety in foggy weather, if we consider *r* = 1 km, *n* = 1.33, and a modulation frequency of *f* =10 MHz, the above rule of thumb indicates that ballistic filtering can be obtained when *µ _{s} ≤* 0.42 m

^{−}^{1}(respectively

*µ*0.084 m

_{s}≤

^{−}^{1}) when

*g*= 0 (respectively

*g*= 0.8). According to the World Meteorological Organization’s recommendations [39], such values of

*µ*correspond to very low visibilities (meteorological optical range

_{s}*MOR*= 3

*/µ*) of 7.5 m (respectively 36 m), still assuming a homogeneous, mostly scattering medium (

_{s}*µ*≫

_{s}*µ*).

#### 5.3. Maximum expectable gain

Lastly, we can estimate the maximum expectable gain under the condition of highly diffusing medium with reduced ballistic contribution (*α* ≪ 1). The gain values derived below may not be quantitatively relevant for practical experiments because many phenomena have been neglected in our analysis so far, such as detector noise, turbulence, spurious ambient illumination, limited dynamics of the detector, etc. Under the above conditions, the maximal expected gain using an intensity modulation scheme and a quadrature demodulation technique is roughly driven by the exponential term, i.e.,
$\mathrm{ln}[{\mathcal{G}}_{bf}]\sim \tau $, which allows one to fairly retrieve the gain values plotted in Fig. 4. As noted above, the value of the maximum expectable gain depends not only on *ω* and *σ* but also on *r* and *µ*.

The evolution of the maximum gain with the physical parameters involved can be analyzed from the expression of *τ*. It is obvious to see that
${\mathcal{G}}_{bf}$ naturally increases with *ω*, but also with *r*. This can be interpreted by noticing that increasing the number of spatial periods along the propagation must increase the efficiency of the ballistic filtering. It is also quite straightforward to show that
${\mathcal{G}}_{bf}$ increases with the absorption coefficient *µ* when the diffusion approximation is valid (*σ/µ* > 10). Though difficult to interpret, it is also possible to demonstrate from the expression of *τ* that
$\mathrm{ln}\left[{\mathcal{G}}_{b{f}_{\mathit{max}}}\right]$ increases with *σ* when *ω/µc* > 8*σ/*3*µ*(1 − *g*)^{2} (which condition is satisfied on the left side of the dotted black lines plotted in Figs. 4(a) and 4(b)), otherwise it should decrease with *σ*.

#### 5.4. Limitations and validity of the study

However, it is important to note here that the diffusion equation for modulated photon transport is admitted to be valid for frequencies lower than *ω/c* < *µ* + *σ* [40], which correspond in Fig. 4 to the region lying below the yellow dashed line. More generally, it is easy to check that when
$g>1-\sqrt{2/3}=0.1835$, the minimum frequency derived above *ω _{min}/c* = 2

*σ/*3(1−

*g*)

^{2}lies outside the domain of validity of the diffusion approximation. This restricts our analysis to situations of quasi-isotropic diffusion regime, in fair agreement with primal premises of the diffusion equation approximation for photon transport.

Finally, we would like to stress again that the quantitative gain values retrieved from Fig. 4 are highly unlikely to be achieved in a practical experiment because we have neglected all sources of experimental imperfections in our work, and, for very high frequencies, the extrapolated gains are obtained from the diffusion approximation beyond its validity. However, the above study makes it possible to understand the interplay of the main physical parameters at stake in ballistic filtering for contrast enhancement.

## 6. Conclusion

We have used the photon diffusion theory and its predictions for transport of intensity modulated light in a diffusing medium along with the noise model for quadrature demodulation scheme, to derive optimal operating points for two application scenarios: scattering parameter estimation as used in diffuse optical imaging and ballistic photon filtering as used in ballistic photon imaging. In the case of estimation of scattering parameters using only diffused light, the Cramer-Rao lower bound on the estimation of scattering parameters was derived and was shown to have a minimum at an optimal modulation frequency. The derived optimal frequency of modulation that achieves minimum variance of estimation depends on the optical penetration depth and the distance of propagation. The evolution of this optimal operating frequency was analytically shown to be increasing with higher absorption coefficient and decreasing with increase in distance of detection and/or the scattering coefficient of the medium. The loss in estimation precision incurred when using a non-optimal operating frequency was quantified. These results pave way for better design of diffuse optical imaging setups that are used in medical instruments. Indeed, given the average scattering parameters of the tissue under test, an optimal frequency can be inferred from the above analytical expressions, and used in the setup to optimize its performance when estimating the parameters of the scattering inhomogeneities.

In the case of ballistic photon imaging using intensity modulated light, an information theoretical metric was introduced to quantify the efficiency of discriminating ballistic photons from diffused photons. A minimum operating angular frequency was derived which appeared to be essentially a linear function of the scattering coefficient of the intervening medium only: it was shown that a significant gain in ballistic filtering can be expected when the angular frequency of modulation *ω* > 2*σc/*3(1 − *g*)^{2}. According to this rule of thumb, ballistic filtering using modulated imaging seems within reach for transport safety applications. On the other hand, the minimum frequencies derived for typical parameters and propagation lengths encountered in biomedical applications seem hardly compatible with available imaging systems, thereby motivating the search for high frequency demodulation imaging systems.

In this approach, diffusion theory of photon transport and noise model of a quadrature de-modulation scheme were tied together using information theoretical tools to provide minimum-variance operating points and to derive the expression of expected gains by taking into account various competing phenomena in the system. It must, however, be kept in mind that the analysis is valid within the applicability of the diffusion approximation, which is supposed to breaks down at large angular modulation frequencies close to *ω/c* = *µ* + *σ*. Moreover, the extremely high numerical values of gain presented in the article are only limited by physical photon noise that is carried forward to the detection scheme. The properties of the detector, like detection noise and detector dynamics were ignored so far to limit the calculations and to focus on the physical limits and optimal operation of such a scheme. However, in real life application, additive detector noise and source phase noise must be taken into account. Incorporating such additional noise factors in our approach is a clear perspective to this work, which will indeed limit the attainable gain values to more realistic values.

## Appendix

## A. Noise model

#### A.1. Noise variance

At each time slice *t _{i}* the photons received at the detector can be modeled as having Poisson noise with variance equal to the mean. Then, the optical noise variance at each slice is

*I*(

*t*) =

_{i}*I*

_{0}(1 +

*M*cos[2

*πt*]). The quadrature components are obtained by weighting each slice with sine and cosine of same frequency. Thus, the total noise variance will propagate as

_{i}/TSimilarly, one shows var $\left(U\right)=\frac{{I}_{0}}{2}$.

#### A.2. Joint PDF of amplitude and phase

Let us consider quadrature components [*U,V*]* ^{T}* as joint Gaussian random variables with mean [

*A*cos[

*ϕ*],

*A*sin[

*ϕ*]]

*and covariance matrix Σ = Diag [Λ*

^{T}^{2},Λ

^{2}]. Then, the joint distribution of the random variables

*U*and

*V*is

For a change of variables, such that
$Z=\sqrt{{U}^{2}+{V}^{2}}$ and Ψ = tan^{−}^{1}[*V/U*] ⇒ *U* = *Z* cos[Ψ], *V* = *Z* sin[Ψ], changed PDF can be obtained by noting that

*u,v*} → {

*z,ψ*} and |.| is the determinant.

## B. Jacobian matrix of the transformation *θ*′ → *θ* for diffuse light only

The Jacobian matrix for *θ*′ → *θ*

## C. Fisher information ${\mathcal{F}}_{\mathcal{D}}(\theta )$

In the presence of diffused light only, the FIM for a change of coordinates *θ*′→*θ* is simply calculated by
${\mathcal{F}}_{\mathcal{D}}(\theta )={J}_{\mathcal{D}}^{T}\mathcal{F}({\theta}^{\prime}){J}_{\mathcal{D}}$

## D. Amplitude detected at ballistic and diffused regions

#### Diffused region $\mathcal{D}$

When the detector (or a pixel) receives only diffused light, the quadrature components detected are given by

#### Ballistic + Diffused region $\mathcal{B}\oplus \mathcal{D}$

When the detector (or a pixel) receives both contributions of diffused light and ballistic light, the quadrature components detected have the following expressions

Similarly,

#### Amplitudes

As a consequence of the above derivations, the amplitudes estimated on detectors that receive only diffused light and detectors that receive both ballistic and diffused light respectively read

## Acknowledgments

This work is funded by CEFIPRA under the project RITFOLD 4604-4 as a collaboration between Institut de Physique de Rennes, Université de Rennes 1, France and Raman Research Institute, Bangalore, India.

We would like to thank Prof. Afshin Daryoush for fruitful discussions. We also acknowledge two anonymous reviewers for their thoughtful suggestions to improve this article.

## References and links

**1. **D. Boas, D. Brooks, E. Miller, C. DiMarzio, M. Kilmer, R. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Proc. Mag. **18**, 57–75 (2001). [CrossRef]

**2. **Y. Y. Schechner and N. Karpel, “Clear underwater vision,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, (CVPR IEEE2004) pp. I-536–I-543 Vol. 1.

**3. **H. Ramachandran and A. Narayanan, “Two-dimensional imaging through turbid media using a continuous wave light source,” Opt. Commun. **154**, 255–260 (1998). [CrossRef]

**4. **W. R. Watkins, D. H. Tofsted, V. G. CuQlock-Knopp, J. B. Jordan, and J. O. Merritt, “Navigation through fog using stereoscopic active imaging,” *Proc. SPIE*4023, Enhanced and Synthetic Vision 2000, 20 (2000). [CrossRef]

**5. **N. Hautiere, J. P. Tarel, and D. Aubert, “Towards Fog-Free In-Vehicle Vision Systems through Contrast Restoration,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, (IEEE, 2007), pp. 1–8.

**6. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**, R41–R93 (1999). [CrossRef]

**7. **M. Brewster and Y. Yamada, “Optical properties of thick, turbid media from picosecond time-resolved light scattering measurements,” Int. J. Heat Mass Tran. **38**, 2569–2581 (1995). [CrossRef]

**8. **B. J. Tromberg, N. Shah, R. Lanning, A. Cerussi, J. Espinoza, T. Pham, L. Svaasand, and J. Butler, “Non-Invasive In Vivo Characterization of Breast Tumors Using Photon Migration Spectroscopy,” Neoplasia **2**, 26–40 (2000). [CrossRef] [PubMed]

**9. **S. B. Colak, M. B. Van DerMark, G. W. Hooft, J. H. Hoogenraad, E. S. Van Der Linden, and F. A. Kuijpers, “Clinical optical tomography and NIR spectroscopy for breast cancer detection,” IEEE J. Sel. Top. Quant. **5**, 1143–1158 (1999). [CrossRef]

**10. **G. Strangman, D. A. Boas, and J. P. Sutton, “Non-invasive neuroimaging using near-infrared light,” Biol. Psychiat. **52**, 679–693 (2002). [CrossRef] [PubMed]

**11. **T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steadystate diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys.19(4) (1992). [CrossRef] [PubMed]

**12. **M. Ghijsen, B. Choi, A. Durkin, S. Gioux, and B. Tromberg, “Real-time simultaneous single snapshot of optical properties and blood flow using coherent spatial frequency domain imaging (cSFDI),” Biomed. Opt. Express **7**, 870–882 (2016). [CrossRef] [PubMed]

**13. **D. A. Benaron and D. K. Stevenson, “Optical time-of-flight and absorbance imaging of biologic media,” Science **259**, 1463–1466 (1993). [CrossRef] [PubMed]

**14. **R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Time-resolved imaging on a realistic tissue phantom: µ(s)’ and µ(a) images versus time-integrated images,” Appl. Opt. **35**, 4533–4540 (1996). [CrossRef] [PubMed]

**15. **E. Gratton, S. Fantini, M. a. Franceschini, G. Gratton, and M. Fabiani, “Measurements of scattering and absorption changes in muscle and brain,” Philos. Trans. R. Soc. B Biol. Sci. **352**, 727–735 (1997). [CrossRef]

**16. **B. J. Tromberg, L. O. Svaasand, T.-T. Tsay, and R. C. Haskell, “Properties of photon density waves in multiple-scattering media,” Appl. Opt. **32**, 607 (1993). [CrossRef] [PubMed]

**17. **H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Simultaneous reconstruction of optical absorption and scattering maps in turbid media from near-infrared frequency-domain data,” Opt. Lett. **20**, 2128–2130 (1995). [CrossRef] [PubMed]

**18. **B. W. Pogue and M. S. Patterson, “Frequency-domain optical absorption spectroscopy of finite tissue volumes using diffusion theory,” Phys. Med. Biol. **39**, 1157–1180 (1994). [CrossRef] [PubMed]

**19. **S. Farsiu, J. Christofferson, B. Eriksson, P. Milanfar, B. Friedlander, A. Shakouri, and R. Nowak, “Statistical detection and imaging of objects hidden in turbid media using ballistic photons,” Appl. Opt. **46**, 5805–5822 (2007). [CrossRef] [PubMed]

**20. **L. Wang, P. P. Ho, C. Liu, G. Zhang, and R. R. Alfano, “Ballistic 2-d imaging through scattering walls using an ultrafast optical kerr gate,” Science **253**, 769–771 (1991). [CrossRef] [PubMed]

**21. **D. Sedarsky, E. Berrocal, and M. Linne, “Quantitative image contrast enhancement in time-gated transillumination of scattering media,” Opt. Express **19**, 1866–1883 (2011). [CrossRef] [PubMed]

**22. **R. Berg, O. Jarlman, and S. Svanberg, “Medical transillumination imaging using short-pulse diode lasers,” Appl. Opt. **32**, 574 (1993). [CrossRef] [PubMed]

**23. **O. Emile, F. Bretenaker, and A. L. Floch, “Rotating polarization imaging in turbid media,” Opt. Lett. **21**, 1706–1708 (1996). [CrossRef] [PubMed]

**24. **L. Mullen, A. Laux, B. Concannon, E. P. Zege, I. L. Katsev, and A. S. Prikhach, “Amplitude-modulated laser imager,” Appl. Optics **43**, 3874–3892 (2004). [CrossRef]

**25. **S. Sudarsanam, J. Mathew, S. Panigrahi, J. Fade, M. Alouini, and H. Ramachandran, “Real-time imaging through strongly scattering media: seeing through turbid media, instantly,” Sci. Rep. **6**, 25033 (2016). [CrossRef] [PubMed]

**26. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28**, 2331–2336 (1989). [CrossRef] [PubMed]

**27. **D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. **36**, 4587 (1997). [CrossRef] [PubMed]

**28. **J. B. Fishkin and E. Gratton, “Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge,” J. Opt. Soc. Am. A. **10**, 127–140 (1993). [CrossRef] [PubMed]

**29. **S. L. Jacques and B. W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt. **13**, 041302 (2008). [CrossRef] [PubMed]

**30. **S. T. Flock, M. S. Patterson, B. C. Wilson, and D. R. Wyman, “Monte Carlo modeling of light propagation in highly scattering tissues. I. Model predictions and comparison with diffusion theory,” IEEE T. Biomed Eng. **36**, 1162–1168 (1989). [CrossRef]

**31. **L. V. Wang and H. Wu, *Biomedical Optics: Principles and Imaging* (Wiley, 2007).

**32. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Wiley, 1999). [CrossRef]

**33. **R. Lange and P. Seitz, “Solid-state time-of-flight range camera,” IEEE J. Quantum Elect. **37**, 390–397 (2001). [CrossRef]

**34. **F. Mufti and R. Mahony, “Statistical analysis of signal measurement in time-of-flight cameras,” ISPRS J. Photogramm. **66**, 720–731 (2011). [CrossRef]

**35. **P. H. Garthwaite, I. T. Jolliffe, and B. Jones, Statistical Inference (Oxford University, 2002).

**36. **H. K. Kim, U. J. Netz, J. Beuthan, and A. H. Hielscher, “Optimal source-modulation frequencies for transport-theory-based optical tomography of small-tissue volumes,” Opt. Express **16**, 18082–18101 (2008). [CrossRef] [PubMed]

**37. **X. Gu, K. Ren, and A. H. Hielscher, “Frequency-domain sensitivity analysis for small imaging domains using the equation of radiative transfer,” Appl. Opt. **46**, 1624–1632 (2007). [CrossRef] [PubMed]

**38. **V. Toronov, E. D’Amico, D. Hueber, E. Gratton, B. Barbieri, and A. Webb, “Optimization of the signal-to-noise ratio of frequency-domain instrumentation for near-infrared spectro-imaging of the human brain,” Opt. Express **11**, 2717–2729 (2003). [CrossRef] [PubMed]

**39. ** WMO-No. 8, “*Guide to Meteorological Instruments and Methods of Observation*,” 7^{th} edition, World Meteo., (2008).

**40. **J. B. Fishkin, S. Fantini, M. J. vande Ven, and E. Gratton, “Gigahertz photon density waves in a turbid medium: Theory and experiments,” Phys. Rev. E **53**, 2307–2319 (1996). [CrossRef]