Theoretical optimal modulation frequencies for scattering parameter estimation and ballistic photon filtering in diffusive media

The efficiency of using intensity modulated light for 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.


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 diffusive photons. A small fraction of the total photons called ballistic photons, undergo forward scattering to reach a detector and they retain the information of its source. It is of wide interest to discriminate the 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]. Diffused 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]. In such cases the time dependent solution to diffusion theory is used to model the transport of pulsed light [11,12] or modulated light [13][14][15][16] 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 temporal frequency-domain photon migration, where an intensity modulated light with modulation frequencies of a few hundred MHz is transmitted through a diffusive 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 [17,18]. 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 forward scattered ballistic photons from multiply scattered (diffuse) photons traversing through a turbid medium. For example, time gated imaging [19,20], polarization gated imaging [21], intensity modulation imaging [22] etc., where information carried by the 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 can also be achieved by an intensity modulation imaging scheme by discriminating them from the diffused light that arrives with reduced modulation index and additional phase.
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 , we provide a brief description of the use of diffusion equation, the transport of modulated light in diffused media and the imaging scheme along with the parameters that will be used in the theoretical analysis. In section , a generalized demodulation scheme is defined and a noise model is presented. Then, in section , the noise model and the diffusion theory are used to calculate the Cramer-Rao bound for parameter estimation. Finally, in section , the ballistic filtering efficiency is analyzed and discussed with respect to real field situations.

IMAGING SCHEME AND DIFFUSION MODEL
FIG. 1: Imaging scheme: A directional source of light with power, P0, forward cone solid angle, Ω, having a limited spectral width (λ ± ∆λ) is detected at a distance r by a detector that subtends an angle dΩ from the source.
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 diffusive medium have been well studied and reviewed [13,[23][24][25]. 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 diffusive media as derived in [14] and many others [25,26] 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.

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 and when detection is carried out sufficiently away from a point source. 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 [27]. Let us now identify some parameters that are important for describing transport in a scattering medium. The scattering length (l σ ) is defined as the mean distance traveled by a photon before the scatter. Its inverse, the reduced scattering coefficient is denoted by σ and defined as σ = µ s (1 − g), where µ s is the scattering coefficient and g, the anisotropy factor, is the mean cosine of the scattering phase function [28,29].
Another important length scale is the transport mean free path (M F P ) denoted by l * , which is the mean distance traveled by photons before they lose their initial directional information. 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 ( µ/D) in diffused medium. Finally, the absorption coefficient (µ) is associated with an absorption length. For ease of reading, the corresponding definitions are tabulated in Tab I.
In addition, for the analysis in the following sections we will use various dimensionless constants that are also given in the right side of Tab I. The dimensionless parameters R δ = r/δ and R * = r/l * correspond to the effective optical attenuation of diffused light and effective optical attenuation of ballistic light, respectively. 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.

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 Effective refractive index n Speed of light in medium c -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, and is given by the Beer's law as I B = ξP 0 e −(µ+σ)r dΩr 2 Ωr 2 = ξP 0 e −(µ+σ)r dΩ Ω . 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 = ξP 0 e −r/δ dΩr 2 4πDr = ξP 0 e −r/δ dΩ 4π r D [14,26]. The expressions of intensities of ballistic and diffused light show that there are clearly two important length scales to be considered, namely, the transport MFP (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 where we use dimensionless variables R * , R δ and Ω = 4π/Ω.

Modulated light in diffused 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 [13,14,24,25]. In this article, we consider an intensity modulated source of light having modulation angular frequency ω, modulation index M and instantaneous intensity i(t) = I 0 (1 + M cos [ωt]). The ballistic light that follows Beer-Lambert's law, is only attenuated and reaches the detector with instantaneous inten- 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 [14]. Then, the instantaneous diffuse light intensity received at the detector is i d (t) = I D (1 + m D cos [ωt + ∆φ]). Without derivation, we present the expression of the reduced modulation index m D and the phase ∆φ which is identically reported in [14,26]: where the parameter q = (1 + 1 + (ω/µc) 2 )/2, is related to the angular frequency of modulation and ranges between [1, ∞) when ω ∈ [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 B /m D 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 II.
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 [14]. 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.

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 = 0 to N ). Then the quadrature components can be estimated by applying sinusoidal weighting of the data points as From the obtained variables U and V , the amplitude A = M I 0 /2 and phase φ of the demodulated signal are respectively given by 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 i in the current notation. Then, a weighted addition of the Poisson distributed data is also Poissonian, indicating that U and V will also be distributed as Poisson random variables with variance I 0 /2 (see Appendix ). For high intensities, the Poisson distribution can be well approximated by a Gaussian distribution (N ) with variance equal to the mean. Thus, we approximate that U and V are distributed as N (A cos [φ] , Λ 2 ) and N (A sin [φ] , Λ 2 ), respectively, with Λ 2 = I 0 /2 and N (x, var(x)) denoting the normal distribution with meanx and variance var(x). Knowing the distribution of the quadrature components, it then possible to derive the joint probability density function (PDF) of the random variables, Z = √ U 2 + V 2 and Ψ = tan −1 (V /U ), by applying the appropriate change of variables to the Gaussian distribution, as given in the Appendix . The expected values of the random variables Z, Ψ and their variance form the parameter vector θ = [A, φ, Λ 2 ], where 2A = M I 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 [30]. Information theoretical studies have been made to accurately determine the phase which contains the depth information of a scene using the above PDF [31]. 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.

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 F(η) ij = − ∂ 2 ln[P (x|η)] ∂ηi∂ηj , 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 just 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.

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, J D , of the transformation θ → θ in the presence of diffuse light only, as denoted by the suffix D. The i th , j th component of this matrix is J D i,j = ∂θ i ∂θ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, ∆φ = R δ 2(q 2 − 1) and Λ 2 D = I D /2 = 3S 0 R * e −R δ /2, where we have set S 0 = ξdΩP 0 /4π. The Jacobian is then calculated to obtain the FIM with respect to parameter vector θ using the transformation F D (θ) = J T D F(θ )J D . Both, the Jacobian and the FIM F D (θ) are shown in Appendices and , 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.

MAXIMAL PRECISION IN ESTIMATION OF SCATTERING PARAMETERS IN DIFFUSE OPTICAL IMAGING
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θ of the parameter vector θ, one has w θθ T w T ≥ w F(θ) −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 a efficient estimator always exists, however, Maximum Likelihood (ML) estimators have been shown to be asymptotically unbiased and efficient for large collection of data, especially for exponential family of distributions [32]. 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 covari-ance bound to reveal some simple conclusions and rules for the estimation problem. We provide the bound as CRB D (θ) = F D (θ) −1 for the invertible matrix F D (θ) which is shown in the Appendix : 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 * increases with its mean value and has an additional frequency dependent term, while the variance of R δ decreases with the R * . The dependence of both the parameters on frequency is the same and has a functional form e (−1+2q)R δ /(−1 + q 2 ). A simple calculus shows that this function reaches a minimum at q opt = 1 + 1 + 4R 2 δ /2R δ . This indicates that there exists an optimal angular frequency ω opt around which the variance of the estimation is minimum. This optimal angular frequency depends only on R δ and is given by the following expression whose evolution is plotted in Fig. 3.a. In the next subsection we analyze the properties and evolution of the optimal frequency ω opt obtained for the estimation of scattering parameters.

Optimal operating frequency
The expression of 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 ω a , is chosen such that ω a /µc = 1 [26]. This is indeed a reasonable choice of operating frequency, justified by the fact that at smaller frequencies on the one hand, the phase change due to diffusion is too small to be detected. On the other hand 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 [26].
To understand the loss incurred by an imperfect choice of operating frequency, we compare the performance of estimation provided by the two frequencies, ω opt and ω a . We plot the ratio of the variance in the estimation of R δ obtained when using ω opt against ω a in Fig. 3.b. The figure provides a quantitative study of the loss in optimal precision in using ω a as opposed to using ω opt . The precision of estimation are equal only when R δ = 5.3, at which point ω opt = ω a , as indicated by the dashed lines in the figures. For example, taking typical values of scattering parameters valid in tissues, such that µ = 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 ω a is independent of the propagation distance 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, higher operating frequencies are suggested as compared to standard operational frequency.
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 µ.  3: (a) Ratio of the optimal frequency ωopt to standard operating frequency ωa = µc as a function of the normalized optical attenuation R δ . (b) Loss in precision in the estimation of R δ using modulation angular frequency ωa as opposed to ωopt as a function of R δ . (c) Contours of optimal frequency of modulation as a function of absorption coefficient µ and reduced scattering coefficient σ and the detection distance r for an intervening medium having refractive index 1. The frequencies can be scaled down by a factor of n for a medium with refractive index n.
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 ω opt = 2µc/ √ R δ = 2 √ µc[r 3(1 + σ/µ)] −1/2 , for sufficiently large value of R δ . 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 [33][34][35].
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 (e.g., the variance of detection 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 spe-cific 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 of ballistic photon filtering in diffusive media, where the goal is to efficiently discriminate the ballistic photons from the diffuse 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.

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 diffusive 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.

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 B⊕D and another set of pixels that receive only diffused light at a region denoted by 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 B ⊕ D as opposed to D. The gain in information, denoted by G bf is thus defined as 11 is the first (upper left) term of the FIM F D (θ) given in Appendix when only diffused light is considered, whereas [F B⊕D (θ)] 11 = [J T B⊕D F(θ )J B⊕D ] 11 stands for the first term of [F B⊕D (θ)] 11 obtained when ballistic and diffused light are simultaneously taken into account. In the above expressions, J B⊕D is the Jacobian of the transformation θ → θ in the presence of ballistic and diffuse light. We do not require here to compute neither the entire Jacobian J B⊕D or the entire FIM F B⊕D (θ) 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 where the A B⊕D and A D are the amplitudes received at the regions B ⊕ D and D, respectively (see Appendix ). Finally, the expression for the gain G bf can be written as where τ = 2(qR δ − R * ), β = e R δ (q−1) and α = I B /I D is the ratio of received ballistic light over diffused light.

Condition for ballistic filtering gain
The contour plot of the ln[G bf ] as a function of the normalized angular modulation frequency (ω/µc) and of σ/µ is shown in Fig. 4, taking Ω = 1 (isotropic emitter) and rµ = 10. 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 coefficient of the medium. The contour for unity gain is not well defined 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 exponential decrease in the amplitude of the cosine term. Moreover, it can be easily checked that it is a sufficient condition to ensure that G bf of Eq. (10a) is greater than unity. The above simple condition can be rewritten as qR δ > R * , which interestingly suggests that 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 simply written as q > R * /R δ = 1 + σ/µ /3, which interestingly appears to be independent of the propagation distance r. The same condition can be expressed in terms of angular frequency ω, yielding the following gain condition This condition is plotted as red dashed line on the contour map 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 ω c > 2 3 σ, and is insensitive to the value of the absorption coefficient µ.
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. The expression derived above can also serve inversely to provide the range of visibilities that can be handled by a ballistic filtering device working at any fixed modulation-demodulation frequency.

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., ln[G bf ] ∼ τ , which allows one to fairly retrieve the gain values plotted in Fig. 4. As noted above, the value of the maximum expectable gain depends on ω and σ but also on r and µ. For example, with µ = 0.002 m −1 , σ/µ = 15 and r = 1 km, the minimum angular frequency is ω/µc = 10. When operating at 1.5 times this frequency one obtains a gain exponent of 5.
The evolution of the maximum gain with the physical parameters involved can be analyzed from the expression of τ . It is obvious to see that 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 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 ln[G bf max ] increases with σ when ω/µc > 8σ/3µ (which condition is plotted in Fig. 4 in dotted black line), otherwise it should decrease with σ. However, it is important to note here that the diffusion equation for modulated photon transport is no more valid for frequencies higher than ω/c > µ + σ [36]. As a result, the minimum frequency condition for achieving ballistic filtering (red dashed line in Fig. 4) remains within the domain of validity, but not the inflection denoted by the black dashed line in Fig. 4.
Once again, we would like to stress 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 obtaind 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. CONCLUSION We have used the photon diffusion theory and its predictions for transport of intensity modulated light in a diffusive 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. The derivation of the optimal operating frequency paves way for better design of diffuse optical imaging se-tups that are used in medical instruments.
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. The theoretical expression of the expected gain at any frequency was also derived and its evolution with respect to physical parameters at hand was discussed.
In this approach, diffusion theory of photon transport and noise model of a quadrature demodulation scheme were tied together using information theoretical tools to provide minimumvariance 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 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.
Jacobian matrix of the transformation θ → θ for diffuse light only The Jacobian matrix for θ → θ Fisher information FD(θ) In the presence of diffused light only, the FIM for a change of coordinates θ → θ is simply calculated by Amplitude detected at ballistic and diffused regions Diffused region D When the detector (or a pixel) receives only diffused light, the quadrature components detected are given by 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