## 1.

## Introduction

An optical inhomogeneity in turbid media (e.g., soft tissue) may imply the existence of a hemorrhage, tumor, or high saturation of oxygenation in a case where absorption and/or scattering characteristics of tissue changes with respect to the concentration of the natural chromophore. For detection and localization of such optical inhomogeneities, illumination with nonionizing electromagnetic (EM) radiation in near-infrared (NIR) or visible region has been utilized. Specifically EM radiation (600 to 1000 nm) is weakly absorbed inside the biological tissue by experiencing many scattering events per unit length, diminishing the availability of detecting ballistic photons.^{1, 2}

NIR radiation is used in diffuse optical imaging (DOI) or tomography (DOT), which is a medical imaging modality to obtain optical parametric maps of living tissue.^{3, 4, 5, 6} The main idea in DOT is the reconstruction of the probed volume by solving a forward model of photon propagation iteratively^{3, 4, 5} and to map spatially resolved spectroscopic chromophore concentrations.^{7, 8, 9, 10} Even though DOT has great advantages, namely use of safe wavelengths, repetitive use, portability, and low cost, it suffers from low spatial resolution with respect to other imaging modalities such as x ray tomography, MRI, and ultrasound imaging.^{11, 12, 13}

Recently, it has been proposed that the location of a defect inside a turbid medium obtained by a direct localization method could be used as *a priori* spatial information for recursive algorithms in DOT.^{14, 15} Direct localization schemes generally depend on transillumination measurements, which can be used for slab, finite-size,
^{14, 16, 17, 18, 19, 20, 21} or circular geometries.^{15} In contrast, for diffuse reflectance measurements, only a few schemes were proposed to improve the accuracy in localization.^{22, 23, 24} Time domain (TD)^{23, 24} and frequency domain (FD) measurements^{22} which are complicated with respect to cw measurements^{11} have attempted to address this issue but no method exists for the fastest and inexpensive method that uses the cw approach.

Knowledge on the depth of an absorber could be incorporated in the process of selecting optimal regularization parameters in both cw and frequency-domain DOI.^{25, 26} Adapting the modulation frequency in accordance to the known depth coordinates allows increased specificity and sensitivity in DOI^{27} and increase in image resolution.^{28, 29, 30, 31, 32} Underdeterminacy of the DOI could also be reduced by introduction of an adaptive grid mesh in the reconstructed volume of interest.^{15, 33} However, co-registration of images from different modalities has subtleness,^{33, 34} and different imaging modalities assess the size of the defects (lesions) differently,^{35} probably due to the inherently different contrast mechanisms in operation.^{33, 36}

The goal of this paper is to recommend a strategy to estimate the depth of an absorber inside a turbid medium. We address the challenge via a single wavelength spatially resolved continuous wave (SRCW) diffuse reflectance for semi-infinite geometry and show that under various probe geometries and optical parameters, it is possible to localize the depth of the absorber to within 0.1-mm precision. The accuracy of the method is analyzed for different combinations of optical properties of the absorber and the medium, as well as for different spatial placement of the cw source and detectors on the surface of the medium.

## 2.

## Theory

An optical defect (absorber) residing in a turbid medium introduces a perturbation to NIR light propagation that is represented as a decrease in measured response. Physically the perturbation is a function of the distance between the defect and the detector when the source and the detector positions are fixed. For a configuration of a single source and multiple detectors in a row (Fig. 1), the perturbation will be less at the farther detectors as long as the defect resides between the source and the closest detector. We propose a method to extract the depth information from the differential effect of perturbation on multiple detectors in a row.

In this study, a solution of time-independent diffusion equation of photon propagation in turbid media is used for simulations.^{37} This solution is derived for the case where the semi-infinite medium is illuminated by a cw source.

The response of a medium containing a spherical inhomogeneity is referred to as the perturbed response [literally total photon-flux-density (*J*
_{T}) in units of detected photons per unit area per unit time at the measurement site]. *J*
_{T} is the summation of two components *J*
_{0} and *J*
_{1}, the response of the unperturbed background medium, and the perturbation introduced by a spherical defect respectively.^{37}

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} J_0 = \delta |\vec E_0 (\vec r)|, \end{equation}\end{document} $${J}_{0}=\delta \left|{\stackrel{\u20d7}{E}}_{0}\left(\stackrel{\u20d7}{r}\right)\right|,$$*x*,

*y*,

*z*) and [TeX:] $\vec E_0 \left({\vec r} \right)$ ${\stackrel{\u20d7}{E}}_{0}\left(\stackrel{\u20d7}{r}\right)$ (in units of photons per unit time per unit volume) is:

^{37}

## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \vec E_0 \left({\vec r} \right) = \frac{{z_0 S_0 }}{{2\pi \delta }}\left({\frac{{3\kappa z\vec r}}{{r^4 }} + \frac{{3z\vec r}}{{r^5 }} {-} \frac{{\kappa \hat z}}{{r^2 }} {-} \frac{{\hat z}}{{r^3 }} {+} \frac{{\kappa ^2 z\vec r}}{{r^3 }}} \right){\rm exp}\left({ - \kappa r} \right){\rm}, \end{equation}\end{document} $${\stackrel{\u20d7}{E}}_{0}\left(\stackrel{\u20d7}{r}\right)=\frac{{z}_{0}{S}_{0}}{2\pi \delta}\left(\frac{3\kappa z\stackrel{\u20d7}{r}}{{r}^{4}}+\frac{3z\stackrel{\u20d7}{r}}{{r}^{5}}-\frac{\kappa \widehat{z}}{{r}^{2}}-\frac{\widehat{z}}{{r}^{3}}+\frac{{\kappa}^{2}z\stackrel{\u20d7}{r}}{{r}^{3}}\right)\mathrm{exp}\left(-\kappa r\right),$$*δ*is the diffusion coefficient which equals [TeX:] $1/({3[ {\mu _\alpha + \mu^ \prime _s } ]})$ $1/\left(3\left[{\mu}_{\alpha}+{\mu}_{s}^{\prime}\right]\right)$ (in unit of length, μ

_{a}and [TeX:] $\mu^ \prime _s$ ${\mu}_{s}^{\prime}$ are absorption and transport scattering coefficients of the medium respectively), κ = (μ

_{a}/δ)

^{1/2}(in unit of length

^{−1}),

*z*

_{0}is the extrapolation depth [TeX:] $({0.7/\mu^ \prime _s })$ $\left(0.7/{\mu}_{s}^{\prime}\right)$ ,

*r*is the magnitude of radius vector [TeX:] $({\vec r})$ $\left(\stackrel{\u20d7}{r}\right)$ , and

*S*

_{0}(in units of photons per unit time) is the photon injection rate which is constant in cw cases.

The perturbation by the spherical defect at the detector site (*J*
_{1}) is:

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \begin{array}{l} \hspace*{-3pt}J_1 = 2\delta q\displaystyle\frac{{\left({\kappa \left| {\vec r_c - \vec r_d } \right| + 1} \right)}}{{\left| {\vec r_c - \vec r_d } \right|^3 }}\exp \left({ - \kappa \left| {\vec r_c - \vec r_d } \right|} \right)\phi _0 \left({\vec r_c } \right) \\[8pt] \hspace*{15pt} - 2\delta p\left\{ {\displaystyle\frac{{[ {\left({\vec r_c - \vec r_d } \right)\,\vec E_0 \left({\vec r} \right)} ]\left({\kappa \left| {\vec r_c - \vec r_d } \right| + 3} \right)z_c }}{{\left| {\vec r_c - \vec r_d } \right|^5 }} - \displaystyle\frac{{E_0 \left({\vec r_c } \right)}}{{\left| {\vec r_c {-} \vec r_d } \right|^3 }}} \right\}\\ \quad\,\,\,\,\,\times\,\exp \left({ - \kappa \left| {\vec r_c {-} \vec r_d } \right|} \right), \end{array}\hspace*{-10pt}\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {J}_{1}=2\delta q\frac{\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|+1\right)}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|}^{3}}\mathrm{exp}\left(-\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|\right){\phi}_{0}\left({\stackrel{\u20d7}{r}}_{c}\right)}\hfill \\ \phantom{\rule{15.0pt}{0ex}}-2\delta p\left\{{\displaystyle \frac{\left[\left({\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right)\phantom{\rule{0.16em}{0ex}}{\stackrel{\u20d7}{E}}_{0}\left(\stackrel{\u20d7}{r}\right)\right]\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|+3\right){z}_{c}}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|}^{5}}-\frac{{E}_{0}\left({\stackrel{\u20d7}{r}}_{c}\right)}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|}^{3}}}\right\}\hfill \\ \phantom{\rule{1em}{0ex}}\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\times \phantom{\rule{0.16em}{0ex}}\mathrm{exp}\left(-\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|\right),\hfill \end{array}\end{array}$$*x*

_{c},

*y*

_{c}, and

*z*

_{c}), [TeX:] $\vec r_d$ ${\stackrel{\u20d7}{r}}_{d}$ is the radius vector of the detector with Cartesian coordinates. The multiplicands

*q*(in units of length) and

*p*(in units of volume) are functions of optical properties of the defect ( [TeX:] $\tilde \mu _a$ ${\stackrel{\u0303}{\mu}}_{a}$ and [TeX:] $\tilde \mu^ \prime _s$ ${\stackrel{\u0303}{\mu}}_{s}^{\prime}$ ), optical properties of the medium ( [TeX:] $\tilde \mu _a$ ${\stackrel{\u0303}{\mu}}_{a}$ and [TeX:] $\tilde \mu^ \prime _s$ ${\stackrel{\u0303}{\mu}}_{s}^{\prime}$ ), and the radius of the defect.

^{37}The

*ϕ*

_{0}( [TeX:] $\vec r_c$ ${\stackrel{\u20d7}{r}}_{c}$ ) (in units of photons per unit area per unit time) is the photon fluence-density-function at the center of the defect.

^{37}

Some of the properties of the multiplicands *q* and *p* are to be stressed. When
[TeX:]
$\tilde \mu _a$
${\stackrel{\u0303}{\mu}}_{a}$
= μ_{a}, *q* is zero. On the other hand, such a relation does not exist for *p* when
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
=
[TeX:]
$\tilde \mu^ \prime _s$
${\stackrel{\u0303}{\mu}}_{s}^{\prime}$
.

It has been reported that an absorber embedded in a semi-infinite medium modifies cw diffuse reflectance such that a shadow on the measurement surface appears.^{38} Methods depending on this idea have been shown to detect the projection of a single defect or two defects onto the measurement surface for the semi-infinite geometry; at the expense of multiple measurements with spatially resolved sampling.

## 3.

## Methods

## 3.1.

### General Approach

## 3.1.1.

#### Generation of ratio-versus-depth curve

As it was reported by Feng, if *q* is nonzero (for
[TeX:]
$\mu _a \ne \tilde \mu _a$
${\mu}_{a}\ne {\stackrel{\u0303}{\mu}}_{a}$
) in Eq. 3, the term with multiplicand *q* is dominant in practice.^{37} When the Eq. 3 is truncated as:

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} J_1 \approx 2\delta q\frac{{\left({\kappa \left| {\vec r_c - \vec r_d } \right| + 1} \right)z_c }}{{\left| {\vec r_c - \vec r_d } \right|}}\exp \left({ - \kappa \left| {\vec r_c - \vec r_d } \right|} \right)\phi _0 \left({\vec r_c } \right) \end{equation}\end{document} $${J}_{1}\approx 2\delta q\frac{\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|+1\right){z}_{c}}{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|}\mathrm{exp}\left(-\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}\right|\right){\phi}_{0}\left({\stackrel{\u20d7}{r}}_{c}\right)$$*J*

_{1}'s calculated for two different detector distances (Fig. 1) is

## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{J_1^P }}{{J_1^D }} = \frac{{2\delta q\frac{{\left({\kappa \left| {\vec r_c - \vec r_d^P } \right| + 1} \right)z_c }}{{\left| {\vec r_c - \vec r_d^P } \right|^3 }}\exp ({ - \kappa | {\vec r_c - \vec r_d^P } |})\phi _0 \left({\vec r_c } \right)}}{{2\delta q\frac{{\left({\kappa \left| {\vec r_c - \vec r_d^D } \right| + 1} \right)z_c }}{{\left| {\vec r_c - \vec r_d^D } \right|^3 }}\exp ({ - \kappa | {\vec r_c - \vec r_d^D } |})\phi \left({\vec r_c } \right)}}, \end{equation}\end{document} $$\frac{{J}_{1}^{P}}{{J}_{1}^{D}}=\frac{2\delta q\frac{\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{P}\right|+1\right){z}_{c}}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{P}\right|}^{3}}\mathrm{exp}\left(-\kappa |{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{P}|\right){\phi}_{0}\left({\stackrel{\u20d7}{r}}_{c}\right)}{2\delta q\frac{\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{D}\right|+1\right){z}_{c}}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{D}\right|}^{3}}\mathrm{exp}\left(-\kappa |{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{D}|\right)\phi \left({\stackrel{\u20d7}{r}}_{c}\right)},$$*P*is for proximal and superscript

*D*is for distal. Since

*q*,

*δ*, and [TeX:] $\phi _0 \left({\vec r_c } \right)$ ${\phi}_{0}\left({\stackrel{\u20d7}{r}}_{c}\right)$ are equal for both the numerator and denominator, they cancel out and Eq. 5 becomes:

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{J_1^P }}{{J_1^D }} = \frac{{\frac{{\left({\kappa \left| {\vec r_c - \vec r_d^P } \right| + 1} \right)}}{{\left| {\vec r_c - \vec r_d^P } \right|^3 }}\exp ({ - \kappa | {\vec r_c - \vec r_d^P } |})}}{{\frac{{\left({\kappa \left| {\vec r_c - \vec r_d^D } \right| + 1} \right)}}{{\left| {\vec r_c - \vec r_d^D } \right|^3 }}\exp ({ - \kappa | {\vec r_c - \vec r_d^D } |})}}. \end{equation}\end{document} $$\frac{{J}_{1}^{P}}{{J}_{1}^{D}}=\frac{\frac{\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{P}\right|+1\right)}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{P}\right|}^{3}}\mathrm{exp}\left(-\kappa |{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{P}|\right)}{\frac{\left(\kappa \left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{D}\right|+1\right)}{{\left|{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{D}\right|}^{3}}\mathrm{exp}\left(-\kappa |{\stackrel{\u20d7}{r}}_{c}-{\stackrel{\u20d7}{r}}_{d}^{D}|\right)}.$$Since the placement of detectors can be arranged, the components of detector vectors (
[TeX:]
$\vec r_d^P$
${\stackrel{\u20d7}{r}}_{d}^{P}$
and
[TeX:]
$\vec r_d^D$
${\stackrel{\u20d7}{r}}_{d}^{D}$
) are fully known. Operands of κ (μ_{a} and
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
) can be assessed by optical methods *in vivo*.
^{39, 40, 41, 42, 43, 44, 45} Please note that *z*
_{c} in Eq. 5 cancels out in Eq. 6. The depth information of the center of the defect remains as a Cartesian component of
[TeX:]
$\vec r_c$
${\stackrel{\u20d7}{r}}_{c}$
with respect to origin – *x*
_{c}, *y*
_{c}, and *z*
_{c}. Assuming that the planar coordinates (*x*
_{c} and *y*
_{c}) of the defect relative to the source are known,^{46} depth (*z*
_{c}) remains the only unknown independent variable on the right-hand side of Eq. 6 implicit in
[TeX:]
$\vec r_c$
${\stackrel{\u20d7}{r}}_{c}$
. A curve generated by Eq. 6 for a range of *z*
_{c} (typically from 5.0 to 40.0 mm) is referred as the ratio-versus-depth curve throughout the text. The generation of ratio-versus-depth curve does not include the knowledge of absorption coefficient of the defect (
[TeX:]
$\tilde \mu _a$
${\stackrel{\u0303}{\mu}}_{a}$
) and radius of the defect (*a*).

In this study, for all cases, the location of proximal detector with respect to source is arranged such that the defect always resides between the source and the proximal detector. Both the source and detectors are taken as extended source and detectors of 1 mm^{2} core-area to have a realistic calculation.

## 3.1.2.

#### Depth estimation

Optically turbid medium is considered to contain an absorbing spherical defect whose scattering coefficient is the same as the medium (Fig. 1). The ratio of perturbed responses of the medium measured at the proximal and distal detectors is:

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} R = \frac{{\left({J_0^P + J_1^P } \right)}}{{\left({J_0^D + J_1^D } \right)}}. \end{equation}\end{document} $$R=\frac{\left({J}_{0}^{P}+{J}_{1}^{P}\right)}{\left({J}_{0}^{D}+{J}_{1}^{D}\right)}.$$*J*

_{0}) can be calculated [Eq. 1] using optical parameters (

*μ*

_{a}and [TeX:] $\mu^{\prime}_s$ ${\mu}_{s}^{\prime}$ ) that are already obtained from look-up tables

^{2}or measured experimentally.

^{39, 40, 41, 42, 43, 44, 45}Then

*J*

_{0}'s of proximal and distal detectors are subtracted in the numerator and denominator of Eq. 8 respectively:

## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} R = \frac{{\left({J_0^P + J_1^P } \right) - J_0^P }}{{\left({J_0^D + J_1^D } \right) - J_0^D }} = \frac{{J_1^P }}{{J_1^D }}. \end{equation}\end{document} $$R=\frac{\left({J}_{0}^{P}+{J}_{1}^{P}\right)-{J}_{0}^{P}}{\left({J}_{0}^{D}+{J}_{1}^{D}\right)-{J}_{0}^{D}}=\frac{{J}_{1}^{P}}{{J}_{1}^{D}}.$$*R*is matched with its corresponding value in the

*y*-axis of ratio-versus-depth curve, the depth is estimated from the

*x*-axis. The ratio-versus-depth curve is generated by Eq. 6 with the same optical parameters used for the calculation of

*J*

_{0}and with the same relative placement of detectors with respect to source and the defect (Fig. 2).

Knowledge on the unperturbed response is required for the application of the proposed method. The method is further extended (described in Sec. 3.4) for the cases where the unperturbed response might be erroneously calculated.

## 3.2.

### Analysis of the Proposed Method

The distance from the source to the planar location of the defect (SC_{x}), source-proximal-detector distance (SPD), and interdetector distance (ID) can be adjusted according to physical conditions (Fig. 1). Absorption coefficient and radius of the defect (*a*) do not take part in the generation of ratio-versus-depth curves [Eq. 6], but still can influence the results since they are present in Eq. 3. Therefore, the analysis of the method comprises two approaches. First, the effect of placement of source and detectors on ratio-versus-depth curves*,* second, the precision of the depth estimation under different probe placement for different defect size and defect/background medium absorption coefficient are analyzed.

The geometry used to derive ratio-versus-depth curves is shown in Fig. 1.

Even if the defect is a pure absorber (
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
=
[TeX:]
$\tilde \mu^ \prime _s$
${\stackrel{\u0303}{\mu}}_{s}^{\prime}$
), the frequency of scattering events inside the defect are increased compared to the background. This could be attributed to light-tissue interaction results from simultaneous scattering and absorption.^{1} The precision of the method is tested quantitatively by calculating error introduced by the approximation [setting *p* = 0 in Eq. 3]. The error is defined as follows: for a particular depth of the center of defect (e.g., point A in Fig. 2), the exact ratio (*y*
_{1}) is calculated by the exact formula [Eq. 3] and plotted as the point (*A*, *y*
_{1}) in Fig. 2. The ratio-versus-depth curve [by Eq. 6] obtained by the same geometry and optical properties is plotted on the same graph. The exact ratio (*y*
_{1}) is projected onto the depth axis to obtain assessed depth via the ratio-versus-depth curve (point B in Fig. 2). Hence, the error at depth *A* is obtained by dividing the absolute value of the difference between *A* and *B* with the diameter of the defect. The divisor is chosen as the diameter of the defect to compare cases with different defect sizes and to avoid depth dependence on the error.

For the quantitative analysis of the method, the errors are calculated for different sets of ID, SC_{x}, *z*
_{c}
_{,} background medium, and defect optical parameters (Fig. 3). Additionally, we studied the error introduced due to possible mistakes in lateral coordinates (Fig. 4). Errors are shown with filled-contour matrices as a function of *z*
_{c} and SC_{x} (Figs. 3, 4).

## Table 1

Parameters used for error calculation in Figs. 3a, 3b, 3c, 3d, 4a, 4b, 4c, 4d. Parameters used for simulations in Fig. 5.

SPD | ID | a (defect radius) | μa (medium) | $\tilde \mu _a$ μ̃a (defect) | $\mu^ \prime _s$ μs′ (medium and defect) | |
---|---|---|---|---|---|---|

(mm) | (mm) | (mm) | (mm−1) | (mm−1) | (mm−1) | |

Figures 3a and 4a | 30.0 | 2.5 | 2.0 | 0.005 | 0.02 | 2.0 |

Figures 3b and 4b | 10.0 | 2.0 | 0.005 | 0.02 | ||

Figures 3c and 4c | 2.5 | 5.0 | 0.005 | 0.02 | ||

Figures 3d, 4d, 5 | 2.5 | 2.0 | 0.050 | 0.20 |

## 3.3.

### Analysis of Probable Errors in Experimental Conditions

The necessity to calculate the response of the background (*J*
_{0}) might be a limiting factor in the applicability of our method for cases where the relatively accurate knowledge on the background optical properties is not available either experimentally or from the literature. A set of optical parameters (predicted from the literature) can be used in iterative algorithms. One approach to use for iterations might be to measure the response of the probed medium with an array of detectors for a configuration in which the defect does not reside between the source and detectors.

The erroneous prediction or estimation of the medium optical parameters is studied and an extension of the method is suggested for improved depth assessment. To analyze this case, the measured response is simulated by feeding a set of optical parameters (actual μ_{a} and
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
) to Eq. 3. Then, the unperturbed response of the medium (*J*
_{0}) is calculated for a range of μ_{a} and
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
selected around the “actual” values and a matrix of calculated ratio values is obtained. The physical constraints require two conditions to be satisfied. First, for each detector (both proximal and distal) the measured response must be less than the corresponding calculated *J*
_{0}. Second, the effect of perturbation must be greater at the proximal detector compared to the distal one (ratio values greater than one). The upper limit of the calculated *J*
_{0} was set not to be higher than three folds of the measured response physically.^{37} The first step is to exclude the set of parameters with a nonphysical outcome from the matrix. The candidate set of optical parameters is then used to estimate the depth of the defect from the corresponding ratio-versus-depth curves for spatially resolved measurements. The next step is to compare the depth values within matrices obtained from different measurements. The depth estimation with smallest standard deviation for the assessed depth values is taken to correspond to the one obtained with the optical values closest to the actual ones (Fig. 5).

## 3.4.

### Experimental Material and Experimental Procedure

An aquarium with darkened sides with dimensions 100.0 × 100.0 × 100.0 mm is filled with 1% intralipid. The aquarium is made of PVC and the top of the aquarium is drilled to hold the 1.1 mm core diameter (whose area is taken as 1.0 mm^{2}) fibers as sources and detectors. A 1-kHz frequency signal is used to operate a 3 mW, 785 nm diode laser (RLD78MA-E Thorlabs). A lock-in amplifier is used to lock the signal to the chosen frequency and eliminate the noise (SRS SR510 Stanford Research). The optical properties of the intralipid solution used in the experiment are μ_{a} = 0.002 mm^{−1} and
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
= 0.827 mm^{−1} by measurement. A 1.0 cm^{3} dark gray cubical rubber is placed at the symmetrical center between source and detector fibers. The source-proximal detector distance is 20.0 mm and the interdetector distance is 10.0 mm to avoid the effect of field-of-view of detectors on measurements.^{47, 48, 49} The source-center of the defect distance is 10.0 mm. Diffuse reflectance values are obtained via descending the rubber by a step size of 1 mm starting at 15 to 40 mm (with respect to the center of the rubber from the sample-ambient layer boundary). The depth range chosen enabled to make measurements of enough intensity. The *J*
_{0} is obtained by measuring spatially-resolved reflectance after the cube is removed from the aquarium.

## 4.

## Results

## 4.1.

### Simulation Results

The diffusely reflected photons injected into a semi-infinite turbid medium by a cw source, traverses the medium in the shape of banana to reach a detector.^{50} The banana-shaped spatial sensitivity (SSP) profile moves toward the medium-air interface with a sharper peak and narrower distribution^{51} when the absorption coefficient of the medium increases. Since the volume of the defect is small with respect to volume probed, the banana-shaped SSP will be used qualitatively in this study.

The optical parameters used in this study are selected to resemble typical *in vivo* parameters for biological tissue cited in literature^{2} (Tables 1, 2). Absorption contrast between the defect and the background medium is fourfold in accordance with previous studies.^{52}

## Table 2

Parameters used for curve generation in Figs. 6, 7, 8.

SCX | SPD | ID | a (defect radius) | μa (medium) | $\tilde \mu _a$ μ̃a (defect) | $\mu^ \prime _s$ μs′ (medium and defect) | |
---|---|---|---|---|---|---|---|

(mm) | (mm) | (mm) | (mm) | (mm−1) | (mm−1) | (mm−1) | |

Fig. 6 | 0.0/5.0 /10.0/15.0/20.0 | 20.0 | 2.5 | ||||

Fig. 7 | 0.0 | 10/20/30 | 2.5 | 3.5 | 0.005 | 0.02 | 2.0 |

Fig. 8 | 0.0 | 20.0 | 2.5/5.0/7.5/10.0 |

The relative placement of detectors and the source with respect to planar position of the defect is an important parameter since different combinations yield different ratio-versus-depth curves (Figs.
6, 7, 8). The effects of different values for source-center of the defect distance (SC_{x}), ID, and SPD are studied in this work (Table 2).

Ratio-versus-depth curves are generated [Eq. 6] for a range of defect depth from 1 to 30 mm with a step size of 1 mm (Figs. 6, 7, 8). To visualize the effect of the neglected term, ratio values obtained by the exact relation given in Eq. 3 are displayed in the same figures (Figs. 6, 7, 8). Used optical and geometrical parameters are displayed in Table 2. For each set of parameter, depths of the center of the defect are taken from 4 to 30 mm with a step size of 2 mm (14 values in Figs. 6, 7, 8).

In a fixed constellation of source relative to the detectors (constant SPD and ID) it is possible to construct ratio-versus-depth curves with different characteristics via changing the placement of the source relative to the planar position of the defect (Fig. 6). As SC_{x} approaches SPD, the perturbation at the proximal detector increases resulting in ratio-versus-depth curves starting at higher ratio values. When a proximal detector is located close to the defect, the scattering effect of the defect is more pronounced and the discrepancy between ratio-versus-depth curves and exact values generated by Eq. 3 become evident for defect depths less than 10 mm. For all cases, as the defect is deeper than 10.0 mm, the discrepancy between ratio-versus-depth curves and exact ratio values diminish due to longer distance to detectors.

The modification introduced in the ratio-versus-depth curve for different values of SPD is given in Fig. 7. The defect is directly below the source for all cases (SC_{x} = 0.0 mm). The ratio-versus-depth curves and exact values overlap quite well for all depths. The discrepancy at depths less than 10.0 mm seems to increase, as SPD gets shorter. A possible explanation would be the augmentation of scattering effect of the defect.

Results for four different values of ID are shown in Fig. 8. As the ID gets longer, the ratio-versus-depth curves starts at higher values because of the decrement of interference of the defect in the detected light intensity of the distal detector. There are minute discrepancies between exact values and ratio-versus-depth curves for all situations that are given.

For the analysis of errors in depth estimation (Fig. 2), the defect is considered to take evenly distributed positions in the region of interest (ROI). ROI is the vertical plane (*y* = 0.0) that spans the distance between the source and the proximal detector in the +*x* direction (0.0 mm ≤ *x* ≤ SPD) and a range for defect depth values in the +*z* direction [ceiling (the radius of the defect) mm + 2.0 mm ≤ *z* ≤ 40.0 mm]. In Figs.
3a, 3b, 3c, 3d, the defect takes 21 evenly spaced locations in the +*x* direction. In Figs.
3a, 3b, 3d the defects are assumed to locate at 37 locations in the +*z* direction spanning from 4.0 to 40.0 mm with a step size of 1.0 mm. For Fig. 3c, however, *z*
_{c} spans from 7.0 to 40.0 mm with 34 depth locations. Optical and geometrical parameters used for error maps are given in Table 1.

In Fig. 3a the error is less than 2.5% of the diameter of the defect when the depth of the center of the defect is below 10.0 mm and its distance from the source is less than 23 mm. The maximum value of error is greater than 15%, which is constrained to a very narrow band at 4.0 mm < *z* < 5.0 mm and 2 mm < *x* < 12 mm.

When ID is extended (from 2.5 to 10.0 mm), the banana of the distal detector becomes more spread [Fig. 3a compared to Fig. 3b]. Therefore, the superficial region of the medium is less equally populated with bananas of proximal and distal detectors resulting in a slightly more extended error in the superficial region. Bigger defect causes [Fig. 3c] errors in the range of 15% when it is in the superficial region (7.0 mm < *z*
_{c} < 8.0 mm and 5 mm < SC_{x} < 23 mm). Increase in the absorption coefficient of the medium shifts bananas of both detectors toward the surface (with a sharper peak and narrower distribution^{51}). Hence, the superficial region becomes more equally populated by the photons of proximal and distal detectors resulting in a decreased error [Fig. 3d compared to Fig. 3a].

In general, for all cases the error is smallest when the source is placed close to the defect (SC_{x} < SPD/2). Photons reaching the two detectors do not equally populate the vertical region below the proximal detector (SC_{x} ∼ SPD), resulting in a higher error in the depth estimation. Increasing ID [Fig. 6b] widens the banana of the distal detector and decreases error. Bigger defects [Fig. 6c] might increase errors when the proximal detector is placed above or close to the defect (error > 0.20), however, the profile of the error distribution is not affected. Error is between 2.5% and 5.0% for the region of *z*
_{c} > 10 mm and SC_{x} < 15 mm.

Errors in the depth estimation due to inaccuracies of lateral positioning of the defect are analyzed for the parameters in Table 1. For error analysis the source is assumed to be directly placed above the center of the defect, however, for each condition the center of the defect is thought to actually be dislocated up to 50% of the radius.

In general, misplacement of the source relative to the lateral position of the defect introduces an error in the depth estimation that is higher for the defects close to the surface and below 20 mm. Placing the source within 0.5 mm from the center of the defect results in lower error (< 5% of the size of the defect) in depth estimation independent of interdetector distances, the size of the defect, or the absorption coefficient of the medium.

The calculated *J*
_{0} is used to obtain the perturbation effect out of the simulated measurements (Fig. 5) (explained in detail in Sec. 3.4).

Results from the extended version of our method suggest that depth estimation of a defect can also be achieved by the use of a range of optical parameters for the calculation of the background response of the probed medium. A relatively large range of optical parameters (0.001 to 0.060 mm^{−1} for *μ*
_{a} and 0.50 to 3.00 mm^{−1} for
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
) is chosen for the simulation that might be limited for real case applications. It should be taken into account that the actual physiologically plausible *μ*
_{a} and
[TeX:]
$\mu^ \prime _s$
${\mu}_{s}^{\prime}$
are taken to reside in the selected range, which would nevertheless need good guesses to be made in advance. Another important fact is to have a finer step size (in the range of 10^{−7}; data not included) for the selected window that would allow converging to the actual optical parameters of the background medium. The necessity is to obtain multiple measurements with various distances from the source (five detectors depicted here) and employ a relatively high computational power (Fig. 5). On the other hand, the ratio-versus-depth curves that necessitate knowledge on the optical properties of the medium are not affected drastically and rough estimation (with an error up to 10%; data not included) would be tolerated in terms of depth estimation.

## 4.2.

### Experimental Results

The comparison of the ratio-versus-depth curves [Eq. 6] and experimental results is shown in Fig. 9. The gray rubber cube of 1.0 cm^{3} is placed in the mid-location of SPD of 20.0 mm (SC_{x} = 10.0 mm) and ID of 10.0 mm. The center of the defect (*z*
_{c}) scanned the depth range of 15.0 to 40.0 mm. For this configuration, the experimental measurements of perturbed response are shown in Fig. 9a in millivolts (mV). Averages of seven independent measurements as a function of depth of the defect and their corresponding standard deviations are presented. The *J*
_{0} measured by the detectors of 20.0 and 30.0 mm distances with respect to source are 0.291 and 0.075 mV respectively. The area of the detectors and the source are taken 1 mm^{2} in ratio-versus-depth curves generation. As expected, the measured perturbed response of proximal detector (distance of 20.0 mm) is higher than that of distal detector [Fig. 9a].

Some of the depths assessed via experimental data and corresponding ratio-versus-depth curves [Fig. 9b] are shown schematically [Fig 9c]. A cross section of the measurement geometry is shown in the figure. The light gray portion is the intralipid solution; the dark gray square shows the cross section of the cubical defect. The black dot inside the dark square is the center-of-mass of the cube, the white one is the assessed depth. Four different depths are represented, 16.0, 24.0, 32.0, and 40.0 mm are chosen.

## 5.

## Discussion

## 5.1.

### General Discussion

In this study, a strategy to optically assess the depth of an absorber inside a semi-infinite turbid medium is proposed. The proposed method involves the use of a single cw source and an array of detectors placed in a row at different distances. By this arrangement a spatially resolved diffuse reflectance is used to obtain the perturbations introduced into the measurements by the defect. The ratio of these perturbations is then checked for the corresponding depth from the ratio-versus-depth curves theoretically constructed with the same optical and spatial parameters. The proposed method does not necessitate knowledge about the size or the absorption coefficient of the defect.

Ratio-versus-depth curves are generated by feeding the optical parameters which can be assessed by either referring to tables^{2} or by experimental means
^{14, 39, 40, 41, 42, 43, 44, 45} and planar coordinates attained by SRCW diffuse reflectance measurements^{38, 46, 52} into Eq. 6. Both experimentally and simulation-wise, the relative placement of the detectors leads to a lower detection at the distal site compared to the proximal one that results in a one-to-one ratio-versus-depth curves approaching one asymptotically.

The desired knowledge on the optical properties of the medium could be a limiting factor for *in vivo* applications.^{53, 54} To minimize this error the unperturbed response of the medium could be obtained experimentally by probing an optically similar region that is expected not to comprise a defect. Additional control schemes are suggested for cases where the optical properties of the background medium are suspected to vary from the look up tables. Our simulations presume that when the calculations for the background response are performed with a range of values around the predicted optical parameters, it is possible to extract the depth of the defect from multiple measurements with an array of detectors.

Notably, the different placements of the source and detectors relative to the planar location of the defect result in curves with different characteristics (e.g., steepness). Eventually changing the arrangement of the detector pair would allow multiple checks for the depth of the defect, enabling gathering of independent values to average for more precise depth estimation. Although the simulations are conducted for a single source-double detector arrangement, the use of detector arrays can increase the precision of the method.

It should be kept in mind that the big defects might influence the error distribution. Generally, the method provides depth estimation with less than 2.5% of the size of the defect (around 0.1 mm precision) for the defects that are deeper than 10 mm from the surface of the medium. This low error is kept independent of the size of the defect [Fig 3c], the absorption coefficient of the defect or of the medium [Fig 3d] as long as the source is placed close to the defect (within 0.5 mm; Fig. 4).

In the experiment, a dark gray cubical rubber is embedded in intralipid solution to simulate a defect in a semi-infinite turbid medium. It was reported by Feng et.al.^{37} that the formulations developed for a spherical defect can be used safely for other types of geometries. Due to low sensitivity of the instruments at our facility, the intralipid solution had to be diluted, resulting in a relatively low absorption coefficient. The simulations on the other hand are performed under more realistic conditions, chosen to resemble the optical properties of biological tissues.^{2} The error analysis also predicts that better results can be obtained in media with higher absorption coefficient compared to the one used for the experiment [Fig. 6d].

Even though it is not reported in detail, when the absorption coefficient of the defect with respect to medium is increased 12-fold instead of 4-fold, the error distribution is similar to those given in this study (Fig. 3). Therefore, the use of a dye can safely be employed to increase the contrast up to 2- or 3-fold.

The experiment to test the proposed method is conducted in an intralipid solution that mimics homogeneous turbid media, which might not be usually the case for biological specimens. However, it can be speculated that the physical principles behind the proposed method should be met when the medium has a degree of heterogeneity, and the background unperturbed response is obtained experimentally and/or validated computationally with iterative methods. For example, the background response of a layered medium can be measured or calculated assuming its optical properties are characterized by a single pair of absorption and reduced scattering coefficients. However, accuracy and limitations of this approach are not studied in this work.

The present method is proposed for a single defect, however, more than one defect could possibly reside in close proximity in the case of biological tissues. When two defects do not coincide vertically, placing the source above one and avoiding the second to reside between source and detectors might allow to assume the second defect as a part of the unperturbed response. Additionally by use of a forward solver of a diffusion equation for two defects, ratio-versus-depth surfaces (two-dimensional version of ratio-versus-depth curves) can be generated to convert the calculated ratio of perturbations into a depth value. On the other hand, for vertically aligned multiple defects the proposed method might be used to assess the average depth but not individual ones. Nevertheless authors agree that the applicability of the method for multiple defects should further be sought.

## 6.

## Conclusions

In this study a cw approach for the localization of an absorber embedded in an otherwise semi-infinite turbid media has been proposed. The motivation behind the study was to develop a method for determination of the depth of the defect independent of the size and absorption coefficient of the defect as much as possible. The method developed here works for semi-infinite geometry, which is a general case in *in vivo* applications. The main contribution of the study is the determination of the depth of the defect with high accuracy (0.1 mm) by using SRCW diffuse reflectance of a single wavelength without iterations when lateral position and background optical properties are predicted. Our method depends on well-known and commonly implemented physical principles^{37} and mathematical derivations. It can be implemented with a NIR device (SRCW, FD, or TD) without any co-registration process. Authors believe that the achieved results could be incorporated in DOI as stated in Sec. 1.

## Acknowledgments

The authors are thankful to the anonymous reviewers for their comments on the earlier version of the manuscript. Burtecin Aksel expresses his deepest thanks to Ayla Aksoy for her comments. This study was supported by the Bogazici University Research Fund through projects 04X102D.