Imaging



Fig. 1
Optical transillumination measurements made with a time-resolved system (left) or an intensity-modulated system (right)




Spectroscopic Measurements



Attenuation of light in the near infrared (NIR) is due to absorption and scattering. The parameter of most interest is absorption which is caused by chromophores of variable concentration such as hemoglobin in its oxygenated and deoxygenated states. In the absence of scattering, the change in light intensity obeys the Beer–Lambert law


$$\displaystyle{ -\ln \frac{I_{{\mathrm{in}}}} {I_{{\mathrm{out}}}} =\mu _{{\mathrm{a}}}\,d =\alpha _{c}[c]d, }$$

(1)
where d is the source-detector separation, which is equal to the optical path length, [c] is the concentration of chromophore c, and α c is the absorption coefficient per unit length per unit concentration of chromophore c and can usually be obtained in vitro. In the presence of scattering, the optical path length of transmitted photons follows a much more complex relationship. Hence, attenuation measurements alone do not allow quantification of chromophore concentration.

Continuous intensity (DC) instruments measure changes in the intensity of light leaving the tissue surface [56]. This is frequently done in a purely spectroscopic manner, i.e., to obtain only global changes in chromophore concentration. In order to quantify concentration changes, additional information is required. One approach is to derive an approximate differential path length factor (DPF) , which restores the approximate Beer–Lambert law for small changes in concentration


$$\displaystyle{ -\delta \ln \frac{I_{{\mathrm{in}}}} {I_{{\mathrm{out}}}} =\mathrm{ DPF\ }\delta \mu _{{\mathrm{a}}}\,d. }$$

(2)
Since there are typically several contributing chromophores, light of different wavelengths in the NIR region is employed, and regression techniques are used to find their relative weightings [26]. It was shown empirically [29] that the DPF is simply the mean time of light multiplied by the speed of light in the tissue. In fact, this relationship follows naturally from the diffusion approximation of light transport [8]. Furthermore, it is equally well approximated by the change in phase of an intensity-modulated system, at least at low modulation frequencies.

Intensity-modulated measurements were first reported by [69]. Most systems use a heterodyne technique to mix the transmitted light with a reference beam of slightly different modulation frequency, thus producing a lower frequency envelope that is easier to detect using RF equipment. Time-resolved systems were first developed using a streak camera [29, 49, 73], an instrument with exceptionally high time resolution in the picosecond range but with high cost, relatively low dynamic range, and a significant inherent temporal nonlinearity due to a sinusoidal ramp voltage. Alternatively, time-correlation single photon counting (TCSPC) systems measure arrival times of individual photons by comparison with a reference pulse using a time-to-amplitude converter (TAC) device [20, 79, 83]. These systems have a high dynamic range and excellent temporal linearity.


Imaging Systems


Imaging methods can be divided into direct systems which seek to detect heterogeneities in tissue by analyzing the transmitted (or, in some cases, reflected) light and indirect systems which attempt to solve the inverse problem of image reconstruction. The latter is the main emphasis of this article although the former is historically the precedent, in a similar manner in which x-ray radiographs were the precursor to x-ray computed tomography (CT).

Transillumination of candle light for a patient suffering from hydrocephalus was reported as early as 1831, but the first significant attempt at diagnostic imaging using optical radiation was for breast lesions and was made by Cutler [27], who used a lamp held under the breast in a darkened room. However, even at this stage, multiple scattering effects caused a notable degradation in image quality. The recognition of this fact led to many attempts to eliminate or minimize the degradation due to scattering ranging from collimation [55] and polarization discrimination [84] to coherence gating using holographic gating [90] or heterodyne detection [82].

With the introduction of time-resolved detectors came the natural attempt to use temporal gating to discriminate early arriving photons (which necessarily have the shortest optical path and therefore suffer the least number of scatterings) from later arriving photons which have undergone multiple scatterings and therefore have ill-determined photon paths. The early implementations of this idea used a Kerr gate as an ultrafast shutter [99]. However, this technique is limited to relatively low-scattering media due to the small dynamic range of the Kerr shutter. Other studies have been based on the streak camera [42] or TCSPC [19] systems described in section “Spectroscopic Measurements.”

The attempt to physically discriminate between photons that have undergone different numbers of scattering events is inherently limited by the statistical likelihood of the low scattering number photons arriving at the detector. For the relatively optically thick tissues that are of interest in breast cancer screening or brain imaging, these photons are overwhelmed by noise. For this reason, indirect methods that solve an inverse problem based on recovering the spatially varying optical parameters that provide the best fit of a photon transport model with the measured data are becoming more prevalent. Within this framework, the three basic strategies (time resolved, intensity modulated, and DC systems) have all been developed and reported. In addition, many different geometrical arrangements have been investigated. Initial studies have been on 2D slice-by-slice imaging, although it is apparent that the photon propagation must in reality be described by a 3D model. Fully 3D methods are now appearing.

In the remainder of this article, we will discuss the inverse problem and the strategies that have been adopted in order to solve it. In order to analyze this problem, we first have to consider the model of photon transport in dense media.



3 Mathematical Modeling and Analysis



Radiative Transfer Equation



In optical imaging, light transport through a medium containing scattering particles is described by transport theory [54]. In transport theory , the particle conservation within a small volume element of phase space is investigated. The wave phenomenon of particles is ignored. The transport theory can be modeled through stochastic methods and deterministic methods. In the stochastic approach, individual particle interactions are modeled as the particles are scattered and absorbed within the medium. The two stochastic methods that have been used in optical imaging are the Monte Carlo method and the random walk theory, of which two, the Monte Carlo is the most often used [5].

In deterministic approach, particle transport is described with integrodifferential equations which can be solved either analytically or numerically [5]. In optical imaging, a widely accepted model for light transport is the radiative transport equation (RTE). The RTE is a one-speed approximation of the transport equation, and thus it basically assumes that the energy (or speed) of the particles does not change in collisions and that the refractive index is constant within the medium. For discussion of photon transport in medium with spatially varying refractive index, see, e.g., [16, 36, 62, 72].

Let $$\Omega \subset \mathbb{R}^{n},\,n = 2\,{\mathrm{or}}\,3$$ denote the physical domain where n is the dimension of the domain. The medium is considered isotropic in the sense that the probability of scattering between two directions depends only on the relative angle between those directions and not on an absolute direction. For discussion of light propagation in anisotropic medium, see, e.g., [44]. Furthermore, let $$\partial \Omega $$ denote the boundary of the domain and $$\hat{s} \in S^{n-1}$$ denote a unit vector in the direction of interest. The RTE is written in time domain as


$$\displaystyle\begin{array}{rcl} & & \frac{1} {c} \frac{\partial \phi (r,\hat{s})} {\partial t} +\hat{ s} \cdot \nabla \phi (r,\hat{s}) + (\mu _{s} +\mu _{a})\phi (r,\hat{s}) \\ & & \quad =\,\mu _{s}\int _{S^{n-1}}\Theta (\hat{s} \cdot \hat{ s}^{{\prime}})\phi (r,\hat{s}^{{\prime}}){\mathrm{d}}\hat{s}^{{\prime}} + q(r,\hat{s}){}\end{array}$$

(3)
and in frequency domain as


$$\displaystyle\begin{array}{rcl} & & \frac{{\mathrm{i}}\omega } {c}\phi (r,\hat{s}) +\hat{ s} \cdot \nabla \phi (r,\hat{s}) + (\mu _{s} +\mu _{a})\phi (r,\hat{s}) \\ & & \quad =\,\mu _{s}\int _{S^{n-1}}\Theta (\hat{s} \cdot \hat{ s}^{{\prime}})\phi (r,\hat{s}^{{\prime}}){\mathrm{d}}\hat{s}^{{\prime}} + q(r,\hat{s}),{}\end{array}$$

(4)
where c is the speed of light in the medium, i is the imaginary unit, ω is the angular modulation frequency of the input signal, and μ s  = μ s (r) and μ a  = μ a (r) are the scattering and absorption coefficients of the medium, respectively. The scattering coefficient represents the probability per unit length of a photon being scattered, and the absorption coefficient represents the probability per unit length of a photon being absorbed. Furthermore, $$\phi (r,\hat{s})$$ is the radiance, $$\Theta (\hat{s} \cdot \hat{ s}^{{\prime}})$$ is the scattering phase function, and $$q(r,\hat{s})$$ is the source inside $$\Omega $$. The radiance can be defined such that the amount of power transfer in the infinitesimal angle $${\mathrm{d}}\hat{s}$$ in direction $$\hat{s}$$ at time t through an infinitesimal area dS is given by


$$\displaystyle{ \phi (r,\hat{s};t)\hat{s} \cdot \hat{\nu }\mathrm{ d}S{\mathrm{d}}\hat{s}, }$$
where $$\hat{\nu }$$ is the normal to the surface dS [54]. The scattering phase function $$\Theta (\hat{s} \cdot \hat{ s}^{{\prime}})$$ describes the probability that a photon with an initial direction $$\hat{s}^{{\prime}}$$ will have a direction $$\hat{s}$$ after a scattering event. In optical imaging, the most usual phase function for isotropic material is the Henyey–Greenstein scattering function [47] which is of the form


$$\displaystyle{ \Theta (\hat{s}\cdot \hat{s}^{{\prime}}) = \left \{\begin{array}{@{}c@{}} \frac{1} {2\pi } \frac{1-g^{2}} {\left (1+g^{2}-2g\hat{s}\cdot \hat{s}^{{\prime}}\right )},\quad n = 2, \\ \frac{1} {4\pi } \frac{1-g^{2}} {\left (1+g^{2}-2g\hat{s}\cdot \hat{s}^{{\prime}}\right )^{3/2}},\quad n = 3, \end{array} \right. }$$

(5)
where g is the scattering shape parameter that defines the shape of the probability density and it gets values between − 1 < g < 1. With the value g = 0, the scattering probability density is a uniform distribution. For forward dominated scattering, g > 0, and for backward dominated scattering, g < 0. The time-domain and frequency-domain representations of the RTE are related through Fourier transform.

In order to obtain a unique solution for the RTE, the ingoing radiance distribution on the boundary $$\partial \Omega $$, that is, $$\phi (r,\hat{s})$$ for $$\hat{s}\cdot \hat{\nu } < 0$$, where $$\hat{\nu }$$ is the outward unit normal, needs to be known [24]. Several boundary conditions can be applied to the RTE [1, 33, 57]. In optical imaging, the boundary condition which assumes that no photons travel in an inward direction at the boundary $$\partial \Omega $$ is used [5]


$$\displaystyle{ \phi (r,\hat{s}) = 0,\quad r \in \partial \Omega,\quad \hat{s} \cdot \hat{n} < 0. }$$

(6)
This boundary condition, also known as the free surface boundary condition and the vacuum boundary condition, implies that once a photon escapes the domain $$\Omega $$, it does not reenter it. The boundary condition (6) can be modified to include a boundary source $$\phi _{0}(r,\hat{s})$$ at the source position $$\varepsilon _{j} \subset \partial \Omega $$, and it can be written in the form [96]


$$\displaystyle{ \phi (r,\hat{s}) = \left \{\begin{array}{@{}c@{\quad }l@{}} \phi _{0}(r,\hat{s}),\quad &r \in \cup _{j}\varepsilon _{j},\quad \hat{s} \cdot \hat{n} < 0 \\ 0, \quad &r \in \partial \Omega \setminus \cup _{j}\varepsilon _{j},\quad \hat{s} \cdot \hat{n} < 0. \end{array} \right. }$$

(7)
In optical imaging, the measurable quantity is the exitance J n (r) on the boundary of the domain. It is defined as [5]


$$\displaystyle{ J_{n}(r) =\int _{S^{n-1}}(\hat{s}\cdot \hat{\nu })\phi (r,\hat{s}){\mathrm{d}}\hat{s},\quad r \in \partial \Omega. }$$

(8)


Diffusion Approximation



In optical imaging, light propagation in tissues is usually modeled with the diffusion approximation (DA) to the RTE. The most typical approach to derive the DA from the RTE is to expand the radiance, the source term, and the phase function into series using the spherical harmonics and truncate the series [5, 24, 33]. If the spherical harmonics series is truncated at the Nth moment, P N approximation is obtained [5, 33]. The first-order spherical harmonics approximation is referred to as the P 1 approximation, and the DA can be regarded as a special case for that. The most typical approach for utilizing the P N approximations in optical imaging has been to use them in angular discretization of the numerical solution of the RTE [14, 100].

An alternative to the P N approximation is the Boltzmann hierarchy approach , in which moments of radiance are used to form a set of coupled equations that approximate the RTE [38]. Furthermore, the DA can be derived using asymptotic techniques [7, 17] leading to generalized diffusion equation or by using projection algebra [33, 57]. If the speed of light is not constant, a diffusion equation with spatially varying indices of refraction can be derived [16].

Here, a short review of the derivation of the DA is given according to [54, 57]. First, the P 1 approximation is derived, and then, the DA is formed as a special case for that. In the DA framework, the radiance is approximated by


$$\displaystyle{ \phi (r,\hat{s}) \approx \frac{1} {\left \vert S^{n-1}\right \vert }\Phi (r) + \frac{n} {\left \vert S^{n-1}\right \vert }\hat{s} \cdot J(r), }$$

(9)
where $$\Phi (r)$$ and J(r) are the photon density and photon current which are defined as


$$\displaystyle\begin{array}{rcl} \Phi (r)& =& \int _{S^{n-1}}\phi (r,\hat{s}){\mathrm{d}}\hat{s}{}\end{array}$$

(10)



$$\displaystyle\begin{array}{rcl} J(r)& =& \int _{S^{n-1}}\hat{s}\phi (r,\hat{s}){\mathrm{d}}\hat{s}.{}\end{array}$$

(11)
By inserting the approximation (9) and similar approximations written for the source term and phase function into Eq. 4 and following the derivation in [5, 54], the P 1 approximation is obtained


$$\displaystyle\begin{array}{rcl} \left ( \frac{{\mathrm{i}}\omega } {c} +\mu _{a}\right )\Phi (r) + \nabla \cdot J(r)& =& q_{0}(r),{}\end{array}$$

(12)



$$\displaystyle\begin{array}{rcl} \left ( \frac{{\mathrm{i}}\omega } {c} +\mu _{a} +\mu _{ s}^{{\prime}}\right )J(r) + \frac{1} {n}\nabla \Phi (r)& =& q_{1}(r),{}\end{array}$$

(13)
where $$\mu _{s}^{{\prime}} = (1 - g_{1})\mu _{s}$$ is the reduced scattering coefficient, q 0(r) and q 1(r) are the isotropic and dipole components of the source, and g 1 is the mean of the cosine of the scattering angle [5, 57]


$$\displaystyle{ g_{1} =\int _{S^{n-1}}(\hat{s} \cdot \hat{s}^{{\prime}})\Theta (\hat{s} \cdot \hat{ s}^{{\prime}}){\mathrm{d}}\hat{s}. }$$

(14)
In the case of the Henyey–Greenstein scattering function , Eq. 5, we have g 1 = g.

To derive the diffusion approximation, it is further assumed that the light source is isotropic, thus q 1(r) = 0, and that $$\frac{{\mathrm{i}}\omega } {c}J(r) = 0$$. The latter assumption, which in time-domain case is of the form $$\frac{1} {c} \frac{\partial J(r)} {\partial t} = 0$$, is usually justified by specifying the condition $$\mu _{a} \ll \mu _{s}^{{\prime}}$$ [5]. Utilizing these approximations, Eq. 13 gives the Fick’s law


$$\displaystyle{ J(r) = -\kappa \nabla \Phi (r), }$$

(15)
where


$$\displaystyle{ \kappa =\kappa (r) = \left (n\left (\mu _{a} +\mu _{ s}^{{\prime}}\right )\right )^{-1} }$$

(16)
is the diffusion coefficient. Substituting Eq. 15 into Eq. 12, the frequency-domain version of the DA is obtained. It is of the form


$$\displaystyle{ -\nabla \cdot \kappa \nabla \Phi (r) +\mu _{a}\Phi (r) + \frac{{\mathrm{i}}\omega } {c}\Phi (r) = q_{0}(r). }$$

(17)
The DA has an analogue in time domain as well. It is of the form


$$\displaystyle{ -\nabla \cdot \kappa \nabla \Phi (r) +\mu _{a}\Phi (r) + \frac{1} {c} \frac{\partial \Phi (r)} {\partial t} = q_{0}(r). }$$

(18)
The time-domain and frequency-domain representations of the DA are related through Fourier transform, similarly as in the case of the RTE.


Boundary Conditions for the DA



The boundary condition (6) cannot be expressed in terms of variables of the diffusion approximation. Instead, there are a few boundary conditions that have been applied to the DA. The simplest boundary condition is the Dirichlet boundary condition which is also referred to as the zero-boundary condition . It sets the photon density to zero on the boundary; thus, $$\Phi (r) = 0,\,r \in \partial \Omega $$ [40, 87]. Alternatively, an extrapolated boundary condition can be used [33, 40, 87]. In the approach, the photon density is set to zero on an extrapolated boundary which is a virtual boundary outside the medium located at a certain distance from the real boundary. Both the zero-boundary condition and the extrapolated boundary condition are physically incorrect, and they have mostly been used because of their mathematical simplicity [40].

The most often used boundary condition in optical imaging is the Robin boundary condition which is also referred to as the partial current boundary condition [4, 24, 33, 40, 54, 87]. It can be derived as follows. Within the P 1 approximation framework (9), the total inward- and outward-directed photon fluxes at a point $$r \in \partial \Omega $$ are


$$\displaystyle\begin{array}{rcl} J^{-}(r)& =& -\int _{\hat{ s}\cdot \hat{n}<0}(\hat{s} \cdot \hat{n})\phi (r,\hat{s}){\mathrm{d}}\hat{s} =\gamma _{n}\Phi (r) -\frac{1} {2}\hat{\nu } \cdot J(r){}\end{array}$$

(19)



$$\displaystyle\begin{array}{rcl} J^{+}(r)& =& \int _{\hat{ s}\cdot \hat{n}>0}(\hat{s} \cdot \hat{n})\phi (r,\hat{s}){\mathrm{d}}\hat{s} =\gamma _{n}\Phi (r) + \frac{1} {2}\hat{\nu } \cdot J(r),{}\end{array}$$” src=”/wp-content/uploads/2016/04/A183156_2_En_21_Chapter_Equ20.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(20)</DIV></DIV>where <SPAN class=EmphasisTypeItalic>γ</SPAN> <SUB><SPAN class=EmphasisTypeItalic>n</SPAN> </SUB>is a dimension-dependent constant which obtains values <SPAN class=EmphasisTypeItalic>γ</SPAN> <SUB>2</SUB> = 1∕<SPAN class=EmphasisTypeItalic>π</SPAN> and <SPAN class=EmphasisTypeItalic>γ</SPAN> <SUB>3</SUB> = 1∕4 [<CITE><A href=57]. To derive the Robin boundary condition for the DA, it is assumed that the total inward-directed photon flux on the boundary is zero; thus,


$$\displaystyle{ J^{-}(r) = 0,\quad r \in \partial \Omega. }$$

(21)
Utilizing Eq. 19 and the Fick’s law (15), the Robin boundary condition can be derived. It is of the form


$$\displaystyle{ \Phi (r) + \frac{1} {2\gamma _{n}}\kappa \frac{\partial \Phi (r)} {\partial \hat{\nu }} = 0,\quad r \in \partial \Omega. }$$

(22)

The boundary condition (22) can be extended to include the reflection on the boundary that is caused by different refractive indices between the object and the surrounding medium. In that case, Eq. 21 is modified to the form


$$\displaystyle{ J^{-}(r) = RJ^{+}(r),\quad r \in \partial \Omega, }$$

(23)
where R = R(x) is the reflection coefficient on the boundary $$\partial \Omega $$, with 0 ≤ R ≤ 1 [57]. Thus, if R = 0, no boundary reflection occurs, and Eq. 23 is reduced into Eq. 21. The parameter R can be derived from Fresnel’s law [40] or, if the refractive index of the surrounding medium is n out = 1, by an experimental fit


$$\displaystyle{ R \approx -1.4399n_{{\mathrm{in}}}^{-2} + 0.7099n_{\mathrm{ in}}^{-1} + 0.6681 + 0.0636n_{\mathrm{ in}}, }$$

(24)
where n in is the refractive index of the medium [87]. Utilizing Eqs. 19 and 20 and the Fick’s law (15), the Robin boundary condition with mismatched refractive indices can be derived. It takes the form


$$\displaystyle{ \Phi (r) + \frac{1} {2\gamma _{n}}\kappa \zeta \frac{\partial \Phi (r)} {\partial \hat{\nu }} = 0,\quad r \in \partial \Omega, }$$

(25)
where ζ = (1 + R)∕(1 − R), with ζ = 1 in the case of no surface reflection. The boundary conditions of the DA for an interface between two highly scattering materials have been discussed, for example, in [4].

The exitance, Eq. 8, can be written utilizing Eqs. 19 and 20, the Fick’s law  (15), and the boundary condition (25). In the DA framework, the exitance is of the form


$$\displaystyle\begin{array}{rcl} J_{n}(r)& =& J^{+}(r) - J^{-}(r) =\hat{\nu } \cdot J(r) \\ & =& -\kappa \frac{\partial \Phi (r)} {\partial \hat{\nu }} = \frac{2\gamma _{n}} {\zeta } \Phi (r),\quad r \in \partial \Omega.{}\end{array}$$

(26)


Source Models for the DA



In the DA framework, light sources are usually modeled by two approximate models, namely, the collimated source model and the diffuse source model. In the case of the collimated source model , the source is modeled as an isotropic point source


$$\displaystyle{ q_{0}(r) =\delta (r - r_{s}), }$$

(27)
where position r s is located at a depth 1∕μ s ′ below the source site [40, 87]. In the case of the diffuse source model, the source is modeled as an inward-directed diffuse boundary current I s at the source position $$\varepsilon _{j} \subset \partial \Omega $$ [87]. In the case of the diffuse source model, Eq. 23 can be modified as


$$\displaystyle{ J^{-}(r) = RJ^{+}(r) + (1 - R)I_{ s},\quad r \in \cup _{j}\varepsilon _{j}. }$$

(28)
Then, following the similar procedure as earlier, the Robin boundary condition with the diffuse source model is obtained. It is of the form


$$\displaystyle{ \Phi (r)+ \frac{1} {2\gamma _{n}}\kappa \zeta \frac{\partial \Phi (r)} {\partial \hat{\nu }} = \left \{\begin{array}{@{}l@{\quad }l@{}} \frac{I_{s}} {\gamma _{n}},\quad &r \in \cup _{j}\varepsilon _{j} \\ 0, \quad &r \in \partial \Omega \setminus \cup _{j}\varepsilon _{j}. \end{array} \right. }$$

(29)


Validity of the DA



The basic condition for the validity of the DA is that the angular distribution of the radiance is almost uniform. In order to achieve this, the medium must be scattering dominated; thus, μ a  ≪ μ s . Most of the tissue types are highly scattering and the DA can be regarded as a good approximation for modeling light propagation within them. The DA has been found to describe light propagation with a good accuracy in situations in which its assumptions are valid [63, 80] and it has been successfully applied in many applications of optical tomography.

However, the condition stating that the angular distribution of the radiance must be almost uniform is violated close to the highly collimated light sources. In addition, the condition cannot be fulfilled in strongly absorbing or low-scattering tissues such as the cerebrospinal fluid which surrounds the brain and fills the brain ventricles. Furthermore, in addition to the above conditions, the DA cannot accommodate realistic boundary conditions or discontinuities at interfaces. The diffusion theory has been found to fail in situations in which its approximations are not valid such as close to the sources [35, 87] and within the low-scattering regions [37, 48].


Numerical Solution Methods for the DA



The analytical solutions of the RTE and its approximations are often restricted to certain specific geometries, and therefore, their exploitability in optical imaging is limited. Therefore, the equations describing light propagation are usually solved with numerical methods. The most often applied numerical methods are the finite difference method and the finite element method (FEM). The latter is generally regarded as more flexible when issues of implementing different boundary conditions and handling complex geometries are considered, and therefore, it is most often chosen as the method for solving equations governing light transport in tissues.

The FE model for the time-varying DA was introduced in [9]. It was later extended to address the topics of boundary conditions and source models [80, 87] and the frequency-domain case of the DA [86]. It can be regarded as the most typical approach to numerically solve the DA.


Hybrid Approaches Utilizing the DA


To overcome the limitations of the diffusion theory close to the light sources and within low-scattering and non-scattering regions, different hybrid approaches and approximate models have been developed.

The hybrid Monte Carlo diffusion method was developed to overcome the limitations of the DA close to the light sources. In the approach, Monte Carlo simulation is combined with the diffusion theory. The method was introduced in [98] to describe light reflectance in a semi-infinite turbid medium, and it was extended for turbid slabs in [97]. In the hybrid Monte Carlo diffusion approach , Monte Carlo method is used to simulate light propagation close to the light source and the DA is analytically solved elsewhere in the domain. Monte Carlo is known to describe light propagation accurately. However, it has the disadvantage of requiring a long computation time. This has effects on computation times of the hybrid Monte Carlo approaches as well. A hybrid radiative transfer–diffusion model to describe light propagation in highly scattering medium was introduced in [93]. In the approach, light propagation is modeled with the RTE close to the light sources, and the DA is used elsewhere in the domain. The solution of the RTE is used to construct a Dirichlet boundary condition for the DA on a fictitious interface within the object. Both the RTE and the DA are numerically solved with the FEM.

Different hybrid approaches and approximate models have been applied for highly scattering media with low-scattering and non-scattering regions. Methods that combine Monte Carlo simulation with diffusion theory have been applied for turbid media with low-scattering regions. The finite element approximation of the DA and the Monte Carlo simulation was combined in [41] to describe light propagation in a scattering medium with a low-scattering layer. However, also in this case, the approach suffers from the time-consuming nature of the Monte Carlo methods. Moreover, the hybrid Monte Carlo diffusion methods often require iterative mapping between the models which increases computation times even more. The radiosity–diffusion model [10, 37] can be applied for highly scattering media with non-scattering regions. The method uses the FE solution of the DA to model light propagation within highly scattering regions and the radiosity model to model light propagation within non-scattering regions. A coupled transport and diffusion model was introduced in [18]. In the model, the transport and diffusion models are coupled, and iterative mapping between the models is used for the forward solution. Furthermore, a coupled radiative transfer equation and diffusion approximation model for optical tomography was introduced in [94] and extended for domains with low-scattering regions in [92]. In the approach, the RTE is used as the forward model in sub-domains in which the assumptions of the DA are not valid and the DA is used elsewhere in the domain. The RTE and DA are coupled through boundary conditions between the RTE and DA sub-domains and solved simultaneously using the FEM.


Green’s Functions and the Robin to Neumann Map


Some insight into light propagation in diffusive media can be gained by examining infinite media. In particular, verification of optical scattering and absorption parameters is frequently made with source and detector fibers immersed in a large container and far from the container walls. In a finite domain, however, we will need to use boundary conditions. We will distinguish between solutions to the homogeneous equation with inhomogeneous boundary conditions  and the inhomogeneous equation with homogeneous boundary conditions. In the latter case, we can use a Green’s function acting on q 0. In the former case, we use a Green’s function acting on a specified boundary function.

We will use the notation $$G_{\Omega }$$ for the Green’s function for the inhomogeneous form of (18) with homogeneous boundary conditions and $$G_{\partial \Omega }$$ for the Green’s function for the homogeneous form of (18) with inhomogeneous boundary conditions, i.e., we have $$G_{\Omega }$$ solving


$$\displaystyle\begin{array}{rcl} & & -\nabla \cdot \kappa (\mathbf{r})\nabla G_{\Omega }(\mathbf{r},\mathbf{r^{{\prime}}},t,t^{{\prime}}) + \left (\mu _{\mathrm{ a}}(\mathbf{r}) + \frac{1} {c} \frac{\partial } {\partial t}\right )G_{\Omega }(\mathbf{r},\mathbf{r^{{\prime}}},t,t^{{\prime}}) =\delta (\mathbf{r^{{\prime}}})\delta (t^{{\prime}}){}\end{array}$$

(30)



$$\displaystyle\begin{array}{rcl} & & \qquad \left.\mathbf{r},\mathbf{r^{{\prime}}}\in \Omega \setminus \partial \Omega,t > t^{{\prime}}\right. \\ & & G_{\Omega }(\mathbf{r_{{\mathrm{d}}}},\mathbf{r^{{\prime}}},t,t^{{\prime}}) + 2\zeta \kappa (\mathbf{r_{\mathrm{ d}}})\frac{\partial G_{\Omega }(\mathbf{r_{{\mathrm{d}}}},\mathbf{r^{{\prime}}},t,t^{{\prime}})} {\partial \nu } = 0 \\ & & \qquad \mathbf{r_{{\mathrm{d}}}} \in \partial \Omega {}\end{array}$$” src=”/wp-content/uploads/2016/04/A183156_2_En_21_Chapter_Equ31.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(31)</DIV></DIV>and <SPAN id=IEq34 class=InlineEquation><IMG alt= solving


$$\displaystyle\begin{array}{rcl} & & -\nabla \cdot \kappa (\mathbf{r})\nabla G_{\partial \Omega }(\mathbf{r},\mathbf{r_{{\mathrm{s}}}},t,t^{{\prime}}) + \left (\mu _{\mathrm{ a}}(\mathbf{r}) + \frac{1} {c} \frac{\partial } {\partial t}\right )G_{\partial \Omega }(\mathbf{r},\mathbf{r_{{\mathrm{s}}}},t,t^{{\prime}}) = 0{}\end{array}$$

(32)



$$\displaystyle\begin{array}{rcl} & & \quad \left.\mathbf{r} \in \Omega \setminus \partial \Omega,t > t^{{\prime}}\right. \\ & & G_{\partial \Omega }(\mathbf{r_{{\mathrm{d}}}},\mathbf{r_{{\mathrm{s}}}},t,t^{{\prime}}) + 2\zeta \kappa (\mathbf{r_{\mathrm{ d}}})\frac{\partial G_{\partial \Omega }(\mathbf{r_{{\mathrm{d}}}},\mathbf{r_{{\mathrm{s}}}},t,t^{{\prime}})} {\partial \nu } =\delta (\mathbf{r_{{\mathrm{s}}}})\delta (t^{{\prime}}) \\ & & \quad \mathbf{r_{{\mathrm{s}}}},\mathbf{r_{{\mathrm{d}}}} \in \partial \Omega. {}\end{array}$$” src=”/wp-content/uploads/2016/04/A183156_2_En_21_Chapter_Equ33.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(33)</DIV></DIV>For a given Green’s function <SPAN class=EmphasisTypeItalic>G</SPAN>, we define the corresponding <SPAN class=EmphasisTypeItalic>Green’s operator</SPAN> as the integral transform with <SPAN class=EmphasisTypeItalic>G</SPAN> as its kernel:<br /> <DIV id=Equb class=Equation><br /> <DIV class=EquationContent><br /> <DIV class=MediaObject><IMG alt=
For the measureable, we define the boundary derivative operator as


$$\displaystyle{\mathcal{B}\,:=\, -\kappa \frac{\partial } {\partial \nu },}$$
where appropriate we will use the simplifying notation


$$\displaystyle{ G^{\mathcal{B}}\,:=\,\mathcal{B}G }$$
to mean the result of taking the boundary data for a Green’s function.

Since (18) is parabolic, we must not simultaneously specify both Dirichlet and Neumann boundary conditions on the whole of $$\partial \Omega $$. The same is true if we convert to the frequency domain and use a complex elliptic equation to describe the propagation of the Fourier transform of $$\Phi $$. Instead, we specify their linear combination through the Robin condition (25). Then, for any specified value q on $$\partial \Omega $$, we will get data y given by (26). The linear mapping $$\Lambda q \rightarrow y$$ is termed the Robin to Neumann map and can be considered the result of a boundary derivative operator $$\mathcal{B}$$ acting on the Green’s operator with kernel $$G_{\partial \Omega }$$


$$\displaystyle{ \Lambda _{{\mathrm{RtN}}}(\kappa,\mu _{{\mathrm{a}}})q = \mathcal{B}\mathcal{G}_{\partial \Omega }q. }$$
Since the Neumann data and Dirichlet data are related by (25), we may also define the Robin to Dirichlet map $$\Lambda _{{\mathrm{RtD}}}(\kappa,\mu _{{\mathrm{a}}})$$ and specify the relationship


$$\displaystyle{ \Lambda _{{\mathrm{RtD}}}(\kappa,\mu _{{\mathrm{a}}}) - 2\zeta \Lambda _{{\mathrm{RtN}}}(\kappa,\mu _{{\mathrm{a}}}) - I = 0 }$$

(34)


The Forward Problem


The Robin to Neumann map is a linear operator mapping boundary sources to boundary data. For the inverse problem, we have to consider a nonlinear mapping from the space of μ a, κ coefficients to the boundary data.

When considering an incoming flux J − with corresponding boundary term q, the data is a function of one variable


$$\displaystyle{ y_{q} = \mathcal{F}_{q}\left (\begin{array}{@{}c@{}} \mu _{{\mathrm{a}}}\\ \kappa \end{array} \right ), }$$

(35)
which gives the boundary data for the particular source term $$q = \mathcal{D}^{-}(J^{-})$$. Using this notation, we consider the forward mapping for a finite number of sources $$\left \{q_{j};j = 1\ldots S\right \}$$ as a parallel set of projections


$$\displaystyle{ \mathbf{y} = \mathbf{\mathcal{F}}\left (\begin{array}{@{}c@{}} \mu _{{\mathrm{a}}}\\ \kappa \end{array} \right ), }$$

(36)
where


$$\displaystyle\begin{array}{rcl} \mathbf{\mathcal{F}}& \,:=\,& \left.(\mathcal{F}_{1},\ldots,\mathcal{F}_{S})\right.^{\mathrm{T}}{}\end{array}$$

(37)



$$\displaystyle\begin{array}{rcl} \mathbf{y}& \,:=\,& \left.(y_{1},\ldots,y_{S})\right.^{\mathrm{T}}.{}\end{array}$$

(38)

We will consider (36) as a mapping from two continuous functions in solution space $$\mu _{{\mathrm{a}}},\kappa \in \mathrm{ X}(\Omega ) \times \mathrm{ X}(\Omega )$$ to continuous functions in data space $$y \in \mathrm{ Y}(\partial \Omega )$$. If the data is sampled as well (which is the case in practice), then $$\mathbf{\mathcal{F}}$$ is sampled at a set of measurement positions $$\left \{\mathbf{r_{{\mathrm{d}}}}_{i};i = 1,\ldots M\right \}$$.

The inverse problem of diffusion-based optical tomography (DOT) is to determine κ, μ a from the values of y for all incoming boundary distributions q. If κ, μ a are found, we can determine μ s ′ through (16).


Schrödinger Form


Problem (17) can be put into Schrödinger form using the Liouville transformation. We make the change of variables $$U =\kappa ^{1/2}\Phi $$, by which (17) becomes


$$\displaystyle\begin{array}{rcl} -\kappa \nabla ^{\,2}\Phi - 2\kappa ^{1/2}\nabla \kappa ^{1/2} \cdot \nabla \Phi + \left (\mu _{\mathrm{ a}} + \frac{{\mathrm{i}}\omega } {c}\right )\Phi = q_{0}.& & {}\\ \end{array}$$
Using


$$\displaystyle{\nabla ^{\,2}U =\kappa ^{1/2}\nabla ^{\,2}\Phi + 2\nabla \Phi \cdot \nabla \kappa ^{1/2} + \Phi \nabla ^{\,2}\kappa ^{1/2}}$$
leads to


$$\displaystyle\begin{array}{rcl} -\nabla ^{\,2}U(\mathbf{r};\omega ) + k^{2}(\mathbf{r};\omega )U(\mathbf{r};\omega ) = \frac{q_{0}(\mathbf{r};\omega )} {\kappa ^{1/2}(\mathbf{r})} & &{}\end{array}$$

(39)



$$\displaystyle\begin{array}{rcl} \quad \mathbf{r} \in \Omega /\partial \Omega & &{}\end{array}$$

(40)



$$\displaystyle\begin{array}{rcl} U(\mathbf{r_{{\mathrm{d}}}};\omega ) + 2\zeta \kappa (\mathbf{r_{{\mathrm{d}}}})\frac{\partial U(\mathbf{r_{{\mathrm{d}}}};\omega )} {\partial \nu } =\kappa ^{1/2}(\mathbf{r_{\mathrm{ d}}})q(\mathbf{r_{{\mathrm{d}}}};\omega )& &{}\end{array}$$

(41)



$$\displaystyle\begin{array}{rcl} \mathbf{r_{{\mathrm{d}}}} \in \partial \Omega,& &{}\end{array}$$

(42)
where


$$\displaystyle{ k^{2} = \frac{\nabla ^{\,2}\kappa ^{1/2}} {\kappa ^{1/2}} + \frac{\mu _{{\mathrm{a}}}} {\kappa } + \frac{{\mathrm{i}}\omega } {c\kappa }. }$$
If k 2 is real (i.e., ω = 0), there exist infinitely many κ, μ a pairs with the same real k 2, so that the measurement of DC data cannot allow the separable unique reconstruction of κ and μ a [6]. For ω ≠ 0, the unique determination of a complex k 2 should be possible by extension of the uniqueness theorem of Sylvester and Uhlmann [91]. From the complex k 2, it is in principle possible to obtain separable reconstruction of first κ from the imaginary part of k 2 and μ a from the real part; see [74] for further discussion.

In a homogeneous medium, with constant optical parameters μ a, κ, we can simplify (42) to


$$\displaystyle\begin{array}{rcl} -\nabla ^{\,2}\Phi (\mathbf{r};\omega ) + k^{2}\Phi (\mathbf{r};\omega ) = q_{ H}(\mathbf{r};\omega ),& &{}\end{array}$$

(43)
with the same boundary condition (25) and with


$$\displaystyle{ k^{2} = \left (\frac{\mu _{{\mathrm{a}}}c +\mathrm{ i}\omega } {c\kappa } \right );\quad q_{H} = \frac{q_{0}(\mathbf{r};\omega )} {\kappa }. }$$

(44)
This equation is also seen directly from (17) for constant κ.

The solution in simple geometries is easily derived using the appropriate Green’s functions [8]. In an infinite medium, this is simply a spherical wave


$$\displaystyle{ \Phi (\mathbf{r};\omega ) \equiv G(\mathbf{r},\mathbf{r_{{\mathrm{s}}}};\omega ) = \frac{{\mathrm{e}}^{\pm k\vert \mathbf{r}-\mathbf{r_{{\mathrm{s}}}}\vert }} {\vert \mathbf{r} -\mathbf{r_{{\mathrm{s}}}}\vert }, }$$

(45)
where the notation G(r, r s; ω) defines Green’s function for a source at position r s. Due to the real part of the wave number ky, this wave is damped. This fact is the main reason that results from diffraction tomography are not always straightforwardly applicable in optical tomography. In particular, for the case ω = 0, the wave is wholly non-propagating. Even as $$\omega \rightarrow \infty $$, the imaginary part of the wave number never exceeds the real part. This is a simple consequence of the parabolic nature of the diffusion approximation. Although hyperbolic approximations can be made too, they do not ameliorate the situation.


Perturbation Analysis



An important tool in scattering problems in general is the approximation of the change in field due to a change in state, developed in a series based on known functions for the reference state. There are two common approaches which we now discuss.


Born Approximation



For the Born approximation, we assume that we have a reference state $${\mathbf{x}}_{0} = \left.(\mu _{{\mathrm{a}}},\kappa )\right.^{\mathrm{T}}$$, with a corresponding wave $$\Phi $$, and that we want to find the scattered wave $$\Phi ^{\delta }$$ due to a change in state $${\mathbf{x}}^{\delta } = \left.(\alpha,\beta )\right.^{\mathrm{T}}$$. We have


$$\displaystyle{ \kappa =\kappa +\beta \,,\quad \mu _{{\mathrm{a}}} =\mu _{{\mathrm{a}}} +\alpha \,. }$$

(46)
Note that it is not necessary to assume that the initial state is homogeneous. Putting (46) into (17) gives


$$\displaystyle\begin{array}{rcl} -\nabla \,\cdot \,(\kappa +\beta )\nabla \tilde{\Phi }(\mathbf{r};\omega ) + \left (\mu _{{\mathrm{a}}} +\alpha + \frac{{\mathrm{i}}\omega } {c}\right )\tilde{\Phi }(\mathbf{r};\omega ) = q_{0}(\mathbf{r};\omega )& &{}\end{array}$$

(47)
with


$$\displaystyle{ \tilde{\Phi } = \Phi + \Phi ^{\delta }. }$$

(48)
Equation 47 can be solved using the Green’s operator for the reference state


$$\displaystyle{ \tilde{\Phi } = \mathcal{G}_{0}\left [q_{0} + \nabla \,\cdot \,\beta \nabla \tilde{\Phi } -\alpha \tilde{ \Phi }\right ]. }$$

(49)
With G 0, Green’s function for the reference state, we have


$$\displaystyle\begin{array}{rcl} \tilde{\Phi }(\mathbf{r};\omega )& =& \Phi (\mathbf{r};\omega ) +\int _{\Omega }\left (G_{0}(\mathbf{r},\mathbf{r^{{\prime}}};\omega )\nabla _{ r^{{\prime}}}\cdot \beta (\mathbf{r^{{\prime}}})\nabla _{ r^{{\prime}}}\tilde{\Phi }(\mathbf{r^{{\prime}}};\omega ) -\alpha (\mathbf{r^{{\prime}}})\tilde{\Phi }(\mathbf{r^{{\prime}}};\omega )\right )\,\mathrm{d^{n}}\mathbf{r}^{{\prime}} \\ & =& \Phi (\mathbf{r};\omega ) -\int _{\Omega }\left (\beta (\mathbf{r^{{\prime}}})\nabla _{ r^{{\prime}}}G_{0}(\mathbf{r},\mathbf{r^{{\prime}}};\omega ) \cdot \nabla _{ r^{{\prime}}}\tilde{\Phi }(\mathbf{r^{{\prime}}};\omega )\right. \\ & & \left.\alpha (\mathbf{r^{{\prime}}})G_{ 0}(\mathbf{r},\mathbf{r^{{\prime}}};\omega )\tilde{\Phi }(\mathbf{r^{{\prime}}};\omega )\right ), {}\end{array}$$

(50)
where we used the divergence theorem and assumed $$\beta (\mathbf{r_{{\mathrm{d}}}}) = 0;\mathbf{r_{{\mathrm{d}}}} \in \partial \Omega $$.

If we define a “potential” as the differential operator


$$\displaystyle{ \mathcal{V}(\alpha,\beta )\,:=\,\nabla \,\cdot \,\beta \nabla -\alpha, }$$

(51)
we can recognize (49) as a Dyson equation and write it in the form


$$\displaystyle{ \left [{\mathsf{I}} -\mathcal{G}_{0}\mathcal{V}\right ]\tilde{\Phi } = \mathcal{G}_{0}q_{0}. }$$

(52)
This may by formally solved by a Neumann series,


$$\displaystyle{ \frac{\mathcal{G}_{0}} {\left [{\mathsf{I}} -\mathcal{G}_{0}\mathcal{V}\right ]} = \mathcal{G}_{0} + \mathcal{G}_{0}\mathcal{V}\mathcal{G}_{0} + \mathcal{G}_{0}\mathcal{V}\mathcal{G}_{0}\mathcal{V}\mathcal{G}_{0} + \cdots }$$

(53)
or, equivalently, by using (48) in (50) to obtain the Born series


$$\displaystyle{ \tilde{\Phi } = \Phi ^{(0)} + \Phi ^{(1)} + \Phi ^{(2)} + \cdots \,, }$$

(54)
where


$$\displaystyle\begin{array}{rcl} \Phi ^{(0)}& =& \Phi {}\\ \Phi ^{(1)}& =& \mathcal{G}_{ 0}\mathcal{V}\Phi {}\\ \Phi ^{(2)}& =& \mathcal{G}_{ 0}\mathcal{V}\mathcal{G}_{0}\mathcal{V}\Phi {}\\ & \vdots & {}\\ \end{array}$$


Rytov Approximation



The Rytov approximation is derived by considering the logarithm of the field as a complex phase [54, 60]:


$$\displaystyle{ \Phi (\mathbf{r};\omega ) =\mathrm{ e}^{u(\mathbf{r};\omega )} }$$

(55)
so that, in place of (48), we have


$$\displaystyle{ \ln \tilde{\Phi } =\ln \Phi + u^{\delta }. }$$

(56)
Putting $$\Phi =\mathrm{ e}^{u_{0}}$$ into (17), we get


$$\displaystyle\begin{array}{rcl} -\Phi \nabla \,\cdot \,\kappa \nabla u_{0} - \Phi \kappa \left \vert \nabla u_{0}\right \vert ^{2} +\tilde{ \Phi }\left (\mu _{\mathrm{ a}} + \frac{{\mathrm{i}}\omega } {c}\right ) = q_{0}.& &{}\end{array}$$

(57)
Putting $$\Phi = \Phi {\mathrm{e}}^{u^{\delta } }$$ and (46) into (17), we get


$$\displaystyle\begin{array}{rcl} & -\tilde{\Phi }\nabla \,\cdot \,(\kappa +\beta )\nabla \left (u_{0} + u^{\delta }\right ) -\tilde{ \Phi }(\kappa +\beta )\left \vert \nabla \left (u_{0} + u^{\delta }\right )\right \vert ^{2} +\tilde{ \Phi }\left (\mu _{{\mathrm{a}}} +\alpha + \frac{{\mathrm{i}}\omega } {c}\right ) = q_{0}.& \\ & & {}\end{array}$$

(58)
Subtracting (57) from (58) and assuming $$\tilde{\Phi } = \Phi $$ over the support of q 0, we get


$$\displaystyle\begin{array}{rcl} & & -\kappa \left (2\nabla u_{0} \cdot \nabla u^{\delta } + \left \vert \nabla u^{\delta }\right \vert ^{2}\right ) -\nabla \,\cdot \,\kappa \nabla u^{\delta } \\ & & = \nabla \,\cdot \,\beta \nabla \left (u_{0} + u^{\delta }\right ) +\beta \left \vert \nabla \left (u_{0} + u^{\delta }\right )\right \vert ^{2} -\alpha.{}\end{array}$$

(59)
We now make use of the relation


$$\displaystyle\begin{array}{rcl} \nabla \,\cdot \,\kappa \nabla u^{\delta }\Phi & =& \Phi \nabla \,\cdot \,\kappa \nabla u^{\delta } + 2\kappa \nabla \Phi \cdot \nabla u^{\delta } + u^{\delta }\nabla \,\cdot \,\kappa \nabla \Phi {}\end{array}$$

(60)



$$\displaystyle\begin{array}{rcl} & =& \Phi \left (\nabla \,\cdot \,\kappa \nabla u^{\delta } + 2\Phi \kappa \nabla u_{0} \cdot \nabla u^{\delta }\right ) + u^{\delta }\nabla \,\cdot \,\kappa \nabla \Phi.{}\end{array}$$

(61)
The last term on the right is substituted from (17) to give


$$\displaystyle\begin{array}{rcl} \nabla \,\cdot \,\kappa \nabla u^{\delta } + 2\kappa \nabla u_{0} \cdot \nabla u^{\delta } = \frac{\nabla \,\cdot \,\kappa \nabla u^{\delta }\Phi } {\Phi } + u^{\delta }\left (\mu _{{\mathrm{a}}} + \frac{{\mathrm{i}}\omega } {c}\right ) -\frac{q_{0}} {\Phi }.& &{}\end{array}$$

(62)
Substituting (62) into (59) and using


$$\displaystyle{\Phi \nabla \,\cdot \,\beta \nabla u_{0} + \Phi \beta \left \vert \nabla u_{0}\right \vert ^{2} = \nabla \,\cdot \,\beta \nabla \left (\Phi u_{ 0}\right )}$$
we arrive at


$$\displaystyle\begin{array}{rcl} & & -\nabla \,\cdot \,\kappa \nabla u^{\delta }\Phi + \left (\mu _{{\mathrm{a}}} + \frac{{\mathrm{i}}\omega } {c}\right )u^{\delta }\Phi \\ & & \quad = \nabla \,\cdot \,\beta \nabla \Phi -\alpha \Phi + \Phi \nabla \,\cdot \,\beta \nabla u^{\delta } +\kappa \left \vert \nabla u^{\delta }\right \vert ^{2}.{}\end{array}$$

(63)
The approximation comes in neglecting the last two terms on the right, which are second order in the small perturbation. The left-hand side is the unperturbed operator, and so the formal solution for $$u^{\delta }\Phi $$ is again obtained through Green’s operator with kernel G 0. Thus, the first-order Rytov approximation becomes


$$\displaystyle\begin{array}{rcl} u^{\delta }(\mathbf{r};\omega )& =& \frac{\Phi ^{(1)}(\mathbf{r};\omega )} {\Phi (\mathbf{r};\omega )}{}\end{array}$$

(64)



$$\displaystyle\begin{array}{rcl} & =& \frac{-1} {\Phi (\mathbf{r};\omega )}\left (\beta (\mathbf{r^{{\prime}}})\nabla _{ r^{{\prime}}}G_{0}(\mathbf{r},\mathbf{r^{{\prime}}};\omega ) \cdot \nabla _{ r^{{\prime}}}\Phi (\mathbf{r^{{\prime}}};\omega )\right. \\ & & \left.\alpha (\mathbf{r^{{\prime}}})G_{ 0}(\mathbf{r},\mathbf{r^{{\prime}}};\omega )\Phi (\mathbf{r^{{\prime}}};\omega )\right ). {}\end{array}$$

(65)

The Rytov approximation is usually argued to be applicable for larger perturbations than the Born approximation, since the neglected terms are small as long as the gradient of the field is slowly varying. See [60] for a much more detailed discussion.

Illustrations of the scattered field in the Born and Rytov formulations are shown in Fig. 2. Since in the frequency domain the field is complex, so is its logarithm. The real part corresponds to the log of the field amplitude and its imaginary part to the phase. From the images in Fig. 2, it is apparent that perturbations are more readily detected in amplitude and phase than in the field itself. This stems from the very high dynamic range of data acquired in optical tomography which in turn stems from the high attenuation and attendant damping. It is the primary motivation for the use of the Rytov approximation, despite the added complications.

A183156_2_En_21_Fig2_HTML.gif


Fig. 2
Top: absorption and scattering images used to generate complex fields. Disk diameter 50 mm, absorption range $$\mu _{{\mathrm{a}}} \in [0.01\mbox{ \textendash }0.04]\,{\mathrm{mm}}^{-1}$$, scatter range $$\mu _{{\mathrm{s}}}^{{\prime}}\in [1\mbox{ \textendash }4]\,{\mathrm{mm}}^{-1}$$. The complex field $$\Phi $$ was calculated using a 2D FEM for a δ-function source on the boundary at the 3 o’clock position. A reference field $$\Phi $$

Only gold members can continue reading. Log In or Register to continue

Related posts:

Segmentation with Shape Priors: Explicit Versus Implicit Representations Transform in Astronomical Data Processing Methods for Multi-dimensional Visual Data Analysis Algorithms from a Non-stochastic Perspective Scattering Impedance Tomography

Stay updated, free articles. Join our Telegram channel
Tags:
Apr 9, 2016 | Posted by in GENERAL RADIOLOGY | Comments Off on Imaging

Full access? Get Clinical Tree
Get Clinical Tree app for offline access