Translator Disclaimer
1 July 2008 Tutorial on diffuse light transport
Abstract
A tutorial introduction to diffuse light transport is presented. The basic analytic equations of time-resolved, steady-state and modulated light transport are introduced. The perturbation method for handling slight heterogeneities in optical properties is outlined. The treatment of boundary conditions such as an air/tissue surface is described. Finite mesh-based numerical methods are introduced to calculate the diffuse light field in complex tissues with arbitrary boundaries. Applications in tissue spectroscopy and imaging illustrate these theoretical and computational tools.

## Introduction

This report is a tutorial introduction to diffuse light transport in biological tissues. Section 1 presents the basics of diffuse light transport, showing the simple equations for time-resolved, steady-state, and modulated light transport. The perturbation method for handling slight heterogeneities in optical properties is introduced. The treatment of an air/tissue surface boundary condition is considered. Section 2 describes numerical methods for simulating light transport in complex tissues. The goal of this work is to provide the novice in biomedical optics with an introduction to diffuse light transport, and the underpinnings of how basic approaches to solving light transport problems can be solved.

## 2.1.

### Optical Properties and Transport Parameters

This tutorial follows the notation of Welch and van Gemert.1 The tissue optical properties are shown in Table 1 . These properties yield the transport properties used in diffusion theory, as given in Table 2 . In discussing light transport, the parameters given in Table 3 are used. The units of power $P$ and light transport $T$ differ for a point source, a line source, and a planar source, but the fluence rate $\varphi =PT$ always has the same units. For a point source, $P$ has units of $W$ and $T$ has units of ${\mathrm{cm}}^{-2}$ . For a line source, $P$ has units of W/cm and $T$ has units of ${\mathrm{cm}}^{-1}$ . For a planar source, $P$ has units of $\mathrm{W}∕{\mathrm{cm}}^{2}$ and $T$ is dimensionless.

## Table 1

Tissue optical properties.

QuantitySymbolUnits
Absorption coefficient ${\mu }_{a}$ ${\mathrm{cm}}^{-1}$
Scattering coefficient ${\mu }_{s}$ ${\mathrm{cm}}^{-1}$
Anisotropy of scattering $g$ dimensionless
Refractive index $n$ dimensionless

## Table 2

Transport properties used in diffusion theory.

QuantitySymbolUnits
Reduced scattering coefficient ${\mu }_{s}^{\prime }={\mu }_{s}\left(1-g\right)$ ${\mathrm{cm}}^{-1}$
Transport mean free path ${\mathit{MFP}}^{\prime }=1∕\left({\mu }_{a}+{\mu }_{s}^{\prime }\right)$ cm
Diffusion length $D={\mathit{MFP}}^{\prime }∕3$ cm
Optical penetration depth $\delta =\sqrt{D∕{\mu }_{a}}$ cm

## Table 3

Parameters for light transport.

QuantitySymbolUnits
Radiant power $P$ W, W/cm, $\mathrm{W}∕{\mathrm{cm}}^{2}$
Radiant energy $Q$ J
Transport $T$ ${\mathrm{cm}}^{-2}$ , ${\mathrm{cm}}^{-1}$ , dimensionless
Fluence rate $\varphi =PT$ $\mathrm{W}∕{\mathrm{cm}}^{2}$
Fluence $\psi =QT$ $\mathrm{J}∕{\mathrm{cm}}^{2}$
Speed of light in tissue $c={c}_{0}∕n$ cm/s
Concentration of radiant energy $C=\varphi ∕c$ $\mathrm{J}∕{\mathrm{cm}}^{3}$

## 2.2.

### Diffuse Versus Nondiffuse Regimes

The movement of light can be discussed in several ways. For example, 1. light moves as an electromagnetic wave in a medium (used in interferometry), 2. light moves as ballistic photons, each with a direction of travel that can be redirected by scattering (used in Monte Carlo simulations), or 3. light can be treated as a concentration of optical energy that diffuses down a concentration gradient (used in diffusion theory simulations). This work is devoted to discussing the third example, the movement of optical energy by diffusion down concentration gradients.

Before proceeding, it is important to understand when diffusion theory is not appropriate. Diffusion assumes that the quantity that is diffusing, in our case optical or radiant energy, does not have a preferential direction of travel. Hence, the net movement of that quantity is by diffusion down a concentration gradient, following Fick’s first law of diffusion (discussed later). When photons are delivered as a collimated beam into a medium, they definitely have a direction of movement. As the photons are scattered by interaction with a tissue, they lose their directionality and hence become eligible for diffusion.

Monte Carlo simulations provide a method for specifying the movement of ballistic photons.2 Figure 1 compares the Monte Carlo and diffusion theory descriptions of the distribution of light in an infinite medium where light is delivered as a collimated beam by an optical fiber submerged in the light-scattering medium. The figure illustrates that in a local region near the source, diffusion theory fails to accurately describe the light distribution. However, distant from the source, diffusion theory is quite good, as if the collimated photons had been thrown forward by the optical fiber and the diffusion process had emanated from this location, a distance $1∕\left[{\mu }_{a}+{\mu }_{s}\left(1-g\right)\right]$ in front of the fiber. This distance is called the transport mean free path, or ${\mathit{MFP}}^{\prime }$ (sometimes called ${l}^{*}$ ). The figure illustrates diffusion theory and the Monte Carlo approach agreement at distances exceeding one ${\mathit{MFP}}^{\prime }$ from the tip of the fiber and exceeding one ${\mathit{MFP}}^{\prime }$ from the apparent source of diffusion located one ${\mathit{MFP}}^{\prime }$ beyond the fiber tip.

## Fig. 1

Monte Carlo simulation of photons launched as collimated beams from optical fibers. The light appears to diffuse from a central point located one ${\mathit{MFP}}^{\prime }=1∕\left[{\mu }_{a}+{\mu }_{s}\left(1-g\right)\right]$ in front of the delivery fiber (red dot at $r=0$ , $z=0.3\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ ). Iso-concentration contours based on Monte Carlo are drawn as black dots and black line. Dashed red circular lines centered around the center of diffusion indicate the prediction of diffusion theory. Near the source, diffusion theory does not agree with the Monte Carlo simulation, but distant from the source agreement is good. The mismatched air/medium boundary at $z=0$ slightly distorts the iso-concentration curves of the Monte Carlo simulation, so Monte Carlo and dashed lines of diffusion theory do not exactly match. (Optical properties: ${\mu }_{a}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu }_{s}=100\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $g=0.90$ , $n=1.4$ .) (Color online only.) The idea of a hybrid Monte Carlo/diffusion theory has been introduced, in which a Monte Carlo simulation launches the photons and diffusion theory propagates them after they have multiply scattered.3 A Monte Carlo simulation launches photons into a scattering medium and propagates the photons according to probability density functions for the step size between photon/tissue interaction sites and the angles of deflection. As the photon propagates initially, the fluence rate in response to a 1-W incident power rate is recorded in an array, ${\varphi }_{\mathit{MC}}\left(x,y,z\right)$ [ $\mathrm{W}∕{\mathrm{cm}}^{2}$ per W], which is the transport ${T}_{\mathit{MC}}\left(x,y,z\right)$ $\left[1∕{\mathrm{cm}}^{2}\right]$ When a photon reaches a depth $z$ that exceeds ${\mathit{MFP}}^{\prime }$ , the Monte Carlo simulation of that photon movement is halted and the photon power is deposited into the local voxel of an array that accumulates the photons as a distributed source $P\left(x,y,z\right)$ [W per voxel]. After launching many photons, the Monte Carlo simulation yields both a ${\varphi }_{\mathit{MC}}\left(x,y,z\right)$ and a distributed source $P\left(x,y,z\right)$ . A point spread function or Green’s function, $G\left({x}^{\prime },{y}^{\prime },{z}^{\prime },x,y,z\right)$ , based on diffusion theory is convolved against $P\left({x}^{\prime },{y}^{\prime },{z}^{\prime }\right)$ to yield the fluence rate, ${\varphi }_{\mathit{DT}}\left(x,y,z\right)=G\otimes P$ , throughout the medium due to diffusion. The final answer is $\varphi ={\varphi }_{\mathit{MC}}+{\varphi }_{\mathit{DT}}$ .

If one ignores the region near a source, then diffusion theory is a very useful and reasonably accurate method for describing the distribution of light in a light-scattering light-absorbing medium. If the observation point is more than one or two transport mean free paths from a source, diffusion theory is usually accurate to within a few percent. Typically in the visible or near-infrared wavelength range, a distance of one ${\mathit{MPF}}^{\prime }$ corresponds to less than 1 or $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . So diffusion theory is quite appropriate and useful when discussing light transport over several millimeters or more.

Diffusion theory also assumes that photons are able to participate in a random walk. Hence, photons should be able to undergo several scattering events before being absorbed. The commonly cited criterion, which is only a rough rule of thumb, is that the ratio ${\mu }_{s}\left(1-g\right)∕{\mu }_{a}$ should exceed 10. If there are sufficient scattering events before an absorption event occurs, then the average photon propagation after many scattering events becomes dependent on ${\mu }_{s}\left(1-g\right)$ rather than on the specific values of ${\mu }_{s}$ and $g$ . The total diffuse reflectance from a tissue that has ${\mu }_{s}\left(1-g\right)∕{\mu }_{a}=10$ is 0.252 (based on Monte Carlo simulations), where the mismatched tissue/air surface boundary has ${n}_{\text{tissue}}∕{n}_{\text{air}}=1.4$ . A study of reflectance versus ${\mu }_{s}\left(1-g\right)∕{\mu }_{a}$ where ${\mu }_{s}$ and $g$ were varied but ${\mu }_{s}\left(1-g\right)$ was conserved, indicated that reflectance becomes independent of $g$ when reflectance exceeds $\sim 0.40$ $\left[{\mu }_{s}\left(1-g\right)∕{\mu }_{a}\phantom{\rule{0.3em}{0ex}}\text{equal}\sim 20\right]$ .4 So the reflectance of a tissue should exceed at least 25% and preferably 40%, or ${\mu }_{s}\left(1-g\right)∕{\mu }_{a}$ should exceed 10 and preferably 20, if one wishes to use diffusion theory reliably. Fortunately, this is usually the case for tissues when using visible to near-infrared wavelengths of light. In tissues with lower reflectance or ratio ${\mu }_{s}\left(1-g\right)∕{\mu }_{a}$ , diffusion theory can still be useful but caution as to the accuracy of diffusion theory is advised.

In summary, diffusion theory is very quick and useful. It is not good near sources. This caution extends to locations near boundaries between regions of differing optical properties, especially the air/tissue surface, and to locations near strongly absorbing objects. Diffusion theory is not so good in tissues with strong absorption, where ${\mu }_{s}\left(1-g\right)∕{\mu }_{a}⪡10$ . For much work in the visible to near-infrared wavelength range, diffusion theory is quite good, and should be understood by every student of biomedical optics.

## 2.3.1.

#### Fluence rate and concentration of optical energy

The fluence rate $\varphi$ $\left[\mathrm{W}∕{\mathrm{cm}}^{2}\right]$ is proportional to the concentration of optical energy $C$ $\left[\mathrm{J}∕{\mathrm{cm}}^{3}\right]$ :

## Eq. 1

$\varphi =cC,$
where $c$ is the speed of light in the medium, $c={c}_{o}∕n$ [cm/s], ${c}_{o}=2.998×{10}^{10}\phantom{\rule{0.3em}{0ex}}\mathrm{cm}∕\mathrm{s}$ , and $n$ is the refractive index of the medium. For example, consider a uniform beam of $1\phantom{\rule{0.3em}{0ex}}\mathrm{W}∕{\mathrm{cm}}^{2}$ irradiance passing through a $1\text{-}{\mathrm{cm}}^{3}$ cube of water (see Fig. 2 ). If the fluence rate $\varphi$ in the nonscattering water $\left(n=1.33\right)$ equals $1\phantom{\rule{0.3em}{0ex}}\mathrm{W}∕{\mathrm{cm}}^{2}$ at $500\text{-}\mathrm{nm}$ wavelength, then the corresponding concentration of photons is $0.18\phantom{\rule{0.3em}{0ex}}\mathrm{pM}$ $\left(0.18×{10}^{-12}\phantom{\rule{0.3em}{0ex}}\text{moles}∕\text{liter}\right)$ :

## Eq. 2

$\varphi \frac{n}{{c}_{0}}\frac{\lambda }{h{c}_{0}}\frac{1000}{{N}_{\mathrm{Av}}}=\left(1\phantom{\rule{0.3em}{0ex}}\mathrm{W}∕{\mathrm{cm}}^{2}\right)\frac{1.33}{{\left(3×{10}^{10}\phantom{\rule{0.3em}{0ex}}\mathrm{cm}∕\mathrm{s}\right)}^{2}}\frac{\left(500×{10}^{-7}\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\right)}{\left(6.63×{10}^{-34}Js\right)}\frac{1000\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{3}∕\text{liter}}{\left(6.02×{10}^{23}\phantom{\rule{0.3em}{0ex}}{\text{mole}}^{-1}\right)}=0.18×{10}^{-12}\mathrm{M}.$
The factor $\lambda ∕hc$ , where $\lambda$ is wavelength, $h$ is Planck’s constant, and $c$ is the in vacuo speed of light, indicates the number of photons per J energy. The factor 1000 indicates the conversion $1000\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{3}∕\text{liter}$ , and ${N}_{\mathrm{Av}}$ is Avagadro’s number for the molecules per mole.

## Fig. 2

The concentration of photons in a cube of water irradiated with $1\phantom{\rule{0.3em}{0ex}}\mathrm{W}∕{\mathrm{cm}}^{2}$ is $0.18\phantom{\rule{0.3em}{0ex}}\mathrm{pM}$ . ## 2.3.2.

Fluence rate can be regarded as a measure of concentration, and the diffusion of light as the movement of photons down concentration gradients. This gradient-driven movement of light is described by Fick’s first law of diffusion, in which the flux $J$ $\left[\mathrm{W}∕{\mathrm{cm}}^{2}\right]$ down a concentration gradient $\partial C∕\partial x$ $\left[\left(\mathrm{J}∕{\mathrm{cm}}^{3}\right)∕\mathrm{cm}\right]$ is

## Eq. 3

$J=-\chi \frac{\partial C}{\partial x}=-D\frac{\partial \varphi }{\partial x},$
where $\chi$ is the diffusivity $\left[{\mathrm{cm}}^{2}∕\mathrm{s}\right]$ equal to $cD$ for light, with the diffusion length $D={\mathit{MFP}}^{\prime }∕3\phantom{\rule{0.3em}{0ex}}\left[\mathrm{cm}\right]$ . In general, $\chi$ can refer to the diffusion of heat, molecules, etc., or radiant energy. For light, $C=\varphi ∕c$ and $\chi =cD$ , and the product $\chi C=\left(cD\right)\left(\varphi ∕c\right)=D\varphi$ . Therefore, the speed of light is canceled in Eq. 3.

The description of the flux of optical energy down a gradient of fluence rate along the dimension $x$ follows Fick’s law, which relates the local gradient to the local flux [Eq. 3]. The description was originally developed to describe neutron scattering in nuclear reactor cores, where the role of absorption was neglected. Consider a small incremental window of cross sectional area $dA$ at $x=0$ , $z=0$ , oriented such that the normal to the window is aligned with the $z$ axis. The flux of light through the window is due to the scatter of light from two small hemispherical regions on either side of the window, which scatter light through the window from both directions. For this example, isotropic scatterers are assumed, described by the reduced scattering coefficient ${\mu }_{s}^{\prime }$ $\left[{\mathrm{cm}}^{-1}\right]$ . The fluence rate upstream from the window along the gradient encounters the scatterers within an incremental volume $dV$ and a power is scattered ${\mu }_{s}^{\prime }\varphi dV$ [W]. A fraction $\mathrm{cos}\left(\theta \right)dA∕\left(4\pi {r}^{2}\right)$ of this scattered power is directed toward the window, which is a distance $r$ from the volume $dV$ , where $\theta$ is the angle between the normal of the window and the direction of scatter toward the window. The fraction of this power aimed at the window that reaches the window is $\mathrm{exp}\left(-r∕{\mathit{MFP}}^{\prime }\right)$ , since scattering and absorption attenuate the light. This attenuation limits the radius of the hemispheres from which the window accepts scattered light. The balance between these two fluxes through the window, the downstream flux minus the upstream flux, yields the net flux through the window. Normalizing the net flux by $dA$ yields the net flux density $J$ $\left[\mathrm{W}∕{\mathrm{cm}}^{2}\right]$ passing through the incremental window:

## Eq. 4

$J={\int }_{z=-\infty }^{0}{\int }_{x=0}^{\infty }{\mu }_{s}^{\prime }\varphi \left(x,z\right)\frac{\mathrm{cos}\left(\theta \right)}{4\pi {r}^{2}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\left(-r∕MF{P}^{\prime }\right)\mathrm{sign}\left(-z\right)2\pi xdxdz.$
The factor $2\pi xdxdz$ is the incremental volume $dV$ expressed as an annular ring. The factor $\mathrm{sign}\left(-z\right)$ causes the integrand to be positive when $z<0$ (upstream of the gradient) and to be negative when $z>0$ (downstream of the gradient). This expression yields a value for the flux density $J$ given a particular gradient $\delta \varphi ∕\delta z$ . The diffusion length is thereby specified as $D=J∕\left(-\delta \varphi ∕\delta z\right)$ , using Fick’s first law.

Fick’s first law holds if the gradient of fluence rate is linear within the two hemispherical regions of integration, i.e., $\varphi \left(z\right)=\varphi \left(0\right)+az$ , where $z$ is the position along the axis perpendicular to the window and $a$ is a constant. This situation is called diffusion theory, and is sometimes called the ${P}^{1}$ approximation. If, however, there is too much curvature in the gradient of $\varphi$ , then the $\varphi$ near the window must be described by higher order terms, $\varphi \left(z\right)=\varphi \left(0\right)+az+b{z}^{2}+\dots$ . An example of this is the so-called ${P}^{3}$ approximation. Such curvature occurs when the point of interest is close to a source of light or a sink for light (an absorber), or near a boundary of mismatched optical properties. Hence, diffusion theory does not work well near discontinuities, and the ${P}^{3}$ approximation is a more appropriate method to describe the transport. However, if the source and point of observation are separated by more than a ${\mathit{MFP}}^{\prime }$ , the ${P}^{1}$ approximation is usually sufficient. The presence of a boundary is discussed later (see Boundaries in see Sec. 2.5).

The use of Fick’s law has two important consequences that should always be kept in mind. First, there is continuity of the fluence rate across any boundary, even if the tissues on either side of the boundary have different optical properties. Second, the flux incident on a boundary equals the flux leaving the boundary, in other words photons do not accumulate on a boundary surface. If two tissues with differing absorption and scattering properties are in contact and light diffuses through their common boundary, the fluence rate $\varphi$ across the boundary will be continuous, although the rate of energy deposition ${\mu }_{a}\varphi$ can be discontinuous.

## 2.3.3.

#### Time-resolved diffusion

Now that photon movement is understood as a diffusion process, consider the standard time-resolved diffusion equation, which is applicable to anything that can diffuse but is now applied to the diffusion of radiant energy. The generic solution to the diffusion equation, for diffusion from a point source, or delta function source, in both space and time is:

## Eq. 5

$C\left(r,t\right)=Q\frac{\mathrm{exp}\left[-{r}^{2}∕\left(4\chi t\right)\right]}{{\left(4\pi \chi t\right)}^{3∕2}}.$
The $Q$ is the point source of whatever substance or energy [units] that will diffuse, for example heat, molecules, photons, etc., which is deposited at the origin $\left(r=0\right)$ at time zero $\left(t=0\right)$ . The second factor has units of ${\mathrm{cm}}^{-3}$ and is referred to as Green’s function. Hence, the concentration $C$ at a distance $r$ from the source at a time $t$ after the initial deposition of $Q$ is expressed as $\left[\text{units}∕{\mathrm{cm}}^{3}\right]$ .

To describe the diffusion of fluence rate $\varphi$ after deposition of a point source of radiant energy, $Q$ $\left[J\right]$ , the speed of light is introduced, since $\varphi =cC$ . Also, the diffusivity $\chi$ is replaced by $cD$ . Finally, the role of absorption is added by including the term $\mathrm{exp}\left(-{\mu }_{a}ct\right)$ , where $ct$ is the path length of individual photons at time $t$ after the impulse, and ${\mu }_{a}$ $\left[{\mathrm{cm}}^{-1}\right]$ is the absorption coefficient. The result is:

## Eq. 6

$\varphi \left(r,t\right)=cQ\frac{\mathrm{exp}\left[-{r}^{2}∕\left(4cDt\right)\right]}{{\left(4\pi cDt\right)}^{3∕2}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\left(-{\mu }_{a}ct\right).$
This is the time-resolved 3-D diffusion equation for light in a medium with scattering and absorption. Figure 3 gives an example of time-resolved diffusion of light from an impulse.

## Fig. 3

Time-resolved diffusion of light from an impulse of light, $Q=1\phantom{\rule{0.3em}{0ex}}\mathrm{J}$ , delivered as a point source at $r=0$ , $t=0$ . The time-resolved fluence rate $\varphi \left(r,t\right)$ is shown at $t=0.1$ , 0.2, and $0.5\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ . ## 2.3.4.

Integration of the prior time-resolved $\varphi \left(r,t\right)$ in Eq. 6 over all time yields the radiant exposure $\psi \left(r\right)$ $\left[\mathrm{J}∕{\mathrm{cm}}^{2}\right]$ , which does not depend on time. The direct integration of Eq. 6 using rules of integration is not immediately obvious. Intuition indicates that $\psi$ should fall exponentially with distance from a source. A numerical integration of $\varphi \left(r,t\right)$ over all time at each position $r$ (not shown here) confirms that at large $r$ , the radiant exposure $\psi \left(r\right)$ falls as $\mathrm{exp}\left(-r∕\text{constant}\right)$ . The constant is given the symbol $\delta$ [cm] and is called the optical penetration depth. The factor $1∕\delta$ is called the effective attenuation coefficient ${\mu }_{\mathrm{eff}}$ . The integration becomes easier by considering conservation of energy. The proper form of $\psi$ is such that the integral of energy deposition over an infinite volume equals the total energy deposited $Q$ . The energy deposition ${\mu }_{a}\psi$ $\left[\mathrm{W}∕{\mathrm{cm}}^{3}\right]$ is integrated over the infinite spatial domain:

## Eq. 7

${\int }_{\infty }{\mu }_{a}\psi \left(r\right)dV={\int }_{r=0}^{\infty }{\mu }_{a}\frac{\mathrm{exp}\left(-r∕\delta \right)}{\text{factor}}4\pi {r}^{2}dr=Q,$
where $dV$ is expressed as an incremental shell volume $4\pi {r}^{2}dr$ . The denominator called factor needs to be specified. The factor needs the term $1∕Q$ to yield the final value $Q$ . The factor needs ${\mu }_{a}$ to cancel the ${\mu }_{a}$ in the numerator. The factor needs $4\pi$ to cancel the $4\pi$ in the numerator. The factor needs $r$ to reduce the ${r}^{2}$ in the numerator to simply $r$ , so that $\mathrm{exp}\left(-r∕\delta \right)rdr$ can be integrated by parts. A factor ${\delta }^{2}$ is generated by this integration, so the factor needs a ${\delta }^{2}$ to cancel this term. In summary, the factor must equal $4\pi {\mu }_{a}{\delta }^{2}r∕Q$ to allow conservation of energy in Eq. 7. Therefore, the expression for the radiant exposure $\psi \left(r\right)$ in response to an impulse of energy $Q$ at $r=0$ and $t=0$ is

## Eq. 8

$\psi \left(r\right)={\int }_{t=0}^{\infty }\varphi \left(r,t\right)dt=Q\frac{\mathrm{exp}\left(-r∕\delta \right)}{4\pi {\mu }_{a}{\delta }^{2}r}=Q\frac{\mathrm{exp}\left(-r∕\delta \right)}{4\pi Dr}.$
The factor ${\mu }_{a}{\delta }^{2}$ is identical to $D$ (derivation not shown here), and the expression using $D$ is also shown in Eq. 8 and is commonly used.

This $\psi \left(r\right)$ is the impulse response of radiant exposure to an impulse of energy deposition $Q$ $\left[J\right]$ . If one has a repetitive sequence of impulses at a frequency $f$ , $\left[1∕\mathrm{s}\right]$ the time averaged power is $P=fQ$ [W]. In the limit of $f\to \text{infinity}$ and $Q\to 0$ such that $P$ is maintained at a constant average power, the superposition of the time-resolved responses $\varphi \left(r,t\right)$ to each impulse $Q$ blurs into a steady-state fluence rate proportional to $\psi \left(r\right)$ , with $Q$ replaced by $P$ :

## Eq. 9

$\varphi \left(r\right)=P\frac{\mathrm{exp}\left(-r∕\delta \right)}{4\pi Dr}.$
This equation describes steady-state diffusion of light from a point source of power $P$ . Figure 4 gives an example of steady-state diffusion of light.

## Fig. 4

Steady-state diffusion of light from a point source $P=1\phantom{\rule{0.3em}{0ex}}\mathrm{W}$ at $r=0$ . The dashed line indicates diffusion with no absorption. The data near $r=0$ are not correct because diffusion theory is inaccurate near the source where the gradient of fluence rate has significant curvature $\left(\mid {\partial }^{2}\varphi ∕\partial {r}^{2}\mid ⪢0\right)$ . (Top) Fluence rate, $\varphi$ . (Bottom) Residual error expressed as $\left(\mathit{MC}\text{-}\mathit{DT}\right)∕DT$ , where $\mathit{MC}$ is the Monte Carlo result and $\mathit{DT}$ is the diffusion theory result. The transport factor $T$ equals $\varphi ∕P$ and characterizes the light transport in a tissue independent of the strength of the source. The following equations list the transport $T$ for 1-D planar diffusion from a planar source ${P}_{\text{planar}}$ $\left[\mathrm{W}∕{\mathrm{cm}}^{2}\right]$ , for 2-D cylindrical diffusion from a line source ${P}_{\mathrm{cyl}}$ [W/cm], and for 3-D spherical diffusion from a point source, ${P}_{\mathrm{sph}}$ [W], derived as outlined before from conservation of energy:

## 10.

${T}_{\text{planar}}\left(z\right)=\frac{\mathrm{exp}\left(-z∕\delta \right)}{2{\mu }_{a}\delta },$
${T}_{\mathrm{cyl}}\left(r\right)\approx \sqrt{\frac{2\delta }{\pi r}}\frac{\mathrm{exp}\left(-r∕\delta \right)}{2\pi {\mu }_{a}{\delta }^{2}},$
${T}_{\mathrm{sph}}\left(r\right)=\frac{\mathrm{exp}\left(-r∕\delta \right)}{4\pi {\mu }_{a}{\delta }^{2}r}.$

## 2.3.5.

#### Modulated diffusion

If the power of a point source is modulated sinusoidally, then a population of photons will be periodically injected into the medium. The time-resolved source is:

## Eq. 11

$P\left(t\right)={P}_{o}\left[1+{M}_{o}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\omega t\right)\right],$
where ${M}_{o}$ is the modulation of the source, $0<{M}_{o}<1$ [dimensionless], and $\omega$ is the angular frequency of modulation [radians/s], $\omega =2\pi f$ , where f is the frequency in [Hz].

These injected photons will diffuse down the gradient from the source toward some point of observation at a detector. The diffusion theory expression for transport from such a time-resolved modulated source (a point source within a medium with no nearby boundaries) is

## 12.

$\varphi \left(r,t\right)={P}_{0}\left[{T}_{ss}\left(r\right)+{M}_{o}\frac{\mathrm{exp}\left(-r{k}^{″}\right)}{4\pi Dr}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\omega t-{k}^{\prime }r\right)\right],$
$={P}_{0}{T}_{ss}\left(r\right)\left[1+{M}_{o}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\left[-r\left({k}^{″}-1∕\delta \right)\right]\mathrm{sin}\left(\omega t-{k}^{\prime }r\right)\right],$
where ${T}_{ss}\left(r\right)$ $\left[1∕{\mathrm{cm}}^{2}\right]$ is the steady-state transport:

## Eq. 13

${T}_{ss}\left(r\right)=\frac{\mathrm{exp}\left(-r∕\delta \right)}{4\pi Dr},$
the factors $k\prime$ and $k″$ are defined:

## 14.

${k}^{″}=\frac{1}{\delta \sqrt{2}}{\left\{{\left[1+{\left(\frac{\omega }{{\mu }_{a}c}\right)}^{2}\right]}^{1∕2}+1\right\}}^{1∕2}=\frac{a}{\delta }\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\left(b\right),$
${k}^{\prime }=\frac{1}{\delta \sqrt{2}}{\left\{{\left[1+{\left(\frac{\omega }{{\mu }_{a}c}\right)}^{2}\right]}^{1∕2}-1\right\}}^{1∕2}=\frac{a}{\delta }\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(b\right),$
where

## 15.

$a={\left[1+{\left(\frac{\omega }{{\mu }_{a}c}\right)}^{2}\right]}^{1∕4},$
$b=\frac{1}{2}\phantom{\rule{0.2em}{0ex}}\mathrm{arctan}\left(\frac{\omega }{{\mu }_{a}c}\right).$
These above expressions are consistent with the appendix of Svaasand 5 The consequence of these expressions is that the fluence rate seen at the detector has a steady-state component ${P}_{o}{T}_{ss}\left(r\right)$ , plus a modulated component characterized by a modulation $M$ [dimensionless] and a phase $\Phi$ [radians]:

## 16.

$\text{Modulation}\phantom{\rule{0.3em}{0ex}}M={M}_{o}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\left[-r\left({k}^{″}-1∕\delta \right)\right],$
$\text{Phase}\phantom{\rule{0.3em}{0ex}}\Phi ={k}^{\prime }r.$
If the angular frequency of modulation $\omega$ is swept over a broad range, then the values of $M$ and $\Phi$ versus frequency constitute a frequency domain description of light transport.

For very low frequencies of modulation as $\omega \to 0$ , the factor $a$ approaches 1 and the factor $b$ approaches 0. Therefore, ${k}^{″}$ approaches $1∕\delta$ and ${k}^{\prime }$ approaches 0. $M$ approaches ${M}_{o}$ , and $\Phi$ approaches 0. So the transport becomes

## Eq. 17

$\varphi \left(r,t\right)={P}_{o}{T}_{ss}\left(r\right)\left[1+{M}_{o}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\omega t\right)\right]=P\left(t\right){T}_{ss}\left(r\right),$
which is simply the steady-state response to the slowly varying light source $P\left(t\right)$ .

For very high frequencies of modulation as $\omega ⪢{\mu }_{a}c$ , the factor $a$ approaches ${\left[\omega ∕\left({\mu }_{a}c\right)\right]}^{1∕2}$ , and the factor $b$ approaches $\pi ∕4$ . Since $\mathrm{cos}\left(\pi ∕4\right)=\mathrm{sin}\left(\pi ∕4\right)=1∕\sqrt{2}$ , the values ${k}^{″}$ and ${k}^{\prime }$ approach the same limiting value:

## Eq. 18

${k}^{″}={k}^{\prime }=\frac{1}{\delta \sqrt{2}}{\left(\frac{\omega }{{\mu }_{a}c}\right)}^{1∕2}.$
In the intermediate regime where $\omega$ is a little below or a little above ${\mu }_{a}c$ , the behavior of $M$ and $\Phi$ varies versus $\omega$ , each in a unique manner, and therefore $M$ and $\Phi$ observed at a particular distance $r$ from the source encode the optical properties ${\mu }_{a}$ and ${\mu }_{s}\left(1-g\right)$ . At very low $\omega ⪡{\mu }_{a}c$ , the phase $\Phi$ is too low to reliably measure. At high $\omega >{\mu }_{a}c$ , the $\Phi$ becomes insensitive to ${\mu }_{a}$ . So the use of modulated light is typically in the intermediate region where $\omega$ is a little less than ${\mu }_{a}c$ . For a tissue with $n=1.37$ and ${\mu }_{a}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , the value $\omega ={\mu }_{a}{c}_{0}∕n$ is $1.65×{10}^{-9}\phantom{\rule{0.3em}{0ex}}\left[\mathrm{rad}∕\mathrm{s}\right]$ , and the frequency $f=\omega ∕\left(2\pi \right)$ is $262\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$ . Typically, modulated diffusion measurements are made with frequencies in the range $f=25\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1000\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$ , which correspond to ${\mu }_{a}=\omega n∕{c}_{0}=0.007\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}0.29\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , respectively. Figure 5 shows $M$ and $\Phi$ versus $r$ for a series of modulation frequencies. Figure 6 shows $M$ and $\Phi$ versus modulation frequency for a series of $r$ positions.

## Fig. 5

Modulation ( $M$ [dimensionless]) and phase ( $\Phi$ [radians]) versus distance ( $r$ [cm]) of observation point from source, for various frequencies ( $f$ [Hz]) of full source modulation, ${M}_{o}=1$ . Optical properties of medium are absorption ${\mu }_{a}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , and reduced scattering ${\mu }_{s}\left(1-g\right)=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ . ## Fig. 6

Modulation $\left(M\right)$ and phase $\left(\Phi \right)$ versus frequency of modulation, observed at $r=1$ , 2, and $3\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ from source. The source modulation and tissue properties are the same as in Fig. 5. ## 2.4.

### Modeling Heterogeneities by the Perturbation Method

Diffusion theory can often provide a quick assessment of the influence of a perturbing object within a tissue. Such a perturbing object can be a region of tissue with different absorption and/or scattering properties, perhaps a tumor with an elevated density of blood vessels due to angiogenesis, a portwine stain lesion with abnormal vasculature, a fibrous growth with increased scattering, or a fluid-filled cyst with decreased scattering.

The perturbation method can use the Green’s function of diffusion theory to assess the influence of a perturbing object on a distant point of observation. The method is introduced using the steady-state diffusion equation, but the method can be developed using frequency-domain transport, as shown in Ref. 6, or time-resolved transport. The method need not be based on diffusion theory, any transport function will suffice.

Consider the example in Fig. 7 . If a power $P$ of light is launched isotropically at a position ${x}_{s},{y}_{s},{z}_{s}$ , the transport to a detector at position ${x}_{d},{y}_{d},{z}_{d}$ over a distance ${r}_{s\to d}$ is

## Eq. 19

${T}_{s\to d}=\frac{\mathrm{exp}\left(-{r}_{s\to d}∕\delta \right)}{4\pi D{r}_{s\to d}},$
and the fluence rate at the detector is ${\varphi }_{d}=P{T}_{s\to d}$ . The background optical properties are ${\mu }_{ao}$ and ${\mu }_{so}^{\prime }$ , yielding a background diffusion constant $D=\left(1∕3\right)∕\left({\mu }_{ao}+{\mu }_{so}^{\prime }\right)$ , and a background optical penetration depth $\delta ={\left(D∕{\mu }_{ao}\right)}^{1∕2}$ .

## Fig. 7

A $0.2\text{-}\mathrm{cm}$ -diam. absorbing object perturbs the field of light from a source, as seen by a detector. The transport of light from the source to the perturbing object ${T}_{s\to p}$ is absorbed by the volume $\left(V\right)$ of the sphere, which has an incremental absorption $\Delta {\mu }_{a}$ above the background tissue optical properties, but the same scattering properties as the background tissue. The power absorbed by the sphere, $\Delta {\mu }_{a}{T}_{s\to p}V$ , has units of power [W], and acts as a negative optical power, ${P}_{p}=-\Delta {\mu }_{a}{T}_{s\to p}V$ , that propagates throughout the tissue and subtracts from the background light supplied by the source. The detector sees a local fluence rate ${\varphi }_{d}=P{T}_{s\to d}+{P}_{p}{T}_{p\to d}$ . Hence, the signal at the detector slightly drops in the presence of the perturbing object. Now consider the presence of a small absorbing object of volume $V$ at position ${x}_{p},{y}_{p},{z}_{p}$ with an extra absorption $\Delta {\mu }_{a}$ , such that its absorption is ${\mu }_{a}={\mu }_{\mathrm{ao}}+\Delta {\mu }_{a}$ . The fluence rate at the object is ${\varphi }_{\mathrm{obj}}=P{T}_{s\to p}$ , where ${T}_{s\to p}$ is computed using a source-object distance ${r}_{s\to p}={\left[{\left({x}_{p}-{x}_{s}\right)}^{2}+{\left({y}_{p}-{y}_{s}\right)}^{2}+{\left({z}_{p}-{z}_{s}\right)}^{2}\right]}^{1∕2}$ . The extra energy deposition in the absorbing object, above the background absorption, is $\Delta {\mu }_{a}{\varphi }_{\mathrm{obj}}$ $\left[\mathrm{W}∕{\mathrm{cm}}^{3}\right]$ . The total extra power deposited in the object is $\Delta {\mu }_{a}{\varphi }_{\mathrm{obj}}V$ [W], which has the units of power. Because the object is absorbing, the object behaves as a source of negative power ${P}_{p}$ :

## Eq. 20

${P}_{p}=-\Delta {\mu }_{a}{\varphi }_{\mathrm{obj}}V.$
This perturbing negative power propagates in a spherically symmetric manner throughout the tissue and reaches the detector, and yields a perturbing fluence rate at the detector, ${P}_{p}{T}_{p\to d}$ . The net fluence rate at the detector is

## Eq. 21

${\varphi }_{d}=P{T}_{s\to d}+{P}_{p}{T}_{p\to d}.$
The first term on the right is positive. The second term on the right is negative for an absorbing object with $\Delta {\mu }_{a}>0$ . If $\Delta {\mu }_{a}<0$ , then ${P}_{p}$ is positive and the perturbation ${P}_{p}{T}_{p\to d}$ is positive.

This is the basic idea of perturbation theory. More details can be found in Ref. 6.

If one has many real sources and many virtual sources due to many perturbing objects, the process is treated as a matrix problem. A vector ${\mathbf{P}}_{j}$ lists all sources. A vector ${\mathbf{P}}_{i}$ lists only the virtual sources, and is a subset of ${\mathbf{P}}_{j}$ . All the ${N}_{0}$ real sources are called ${\mathbf{P}}_{j}\left(j=1:{N}_{0}\right)$ , and the ${N}_{i}$ virtual sources due to the perturbing objects are called ${\mathbf{P}}_{i}\left(i=1:{N}_{i}\right)={\mathbf{P}}_{j}\left(j={N}_{0}+1:{N}_{j}\right)$ . A perturbation matrix $\mathbf{M}$ is defined with ${N}_{j}$ columns to match ${\mathbf{P}}_{j}$ and ${N}_{i}$ rows to match ${\mathbf{P}}_{i}$ . The $\mathbf{M}$ matrix consists of elements that scale each source in ${\mathbf{P}}_{j}$ to yield a contribution to an updated value of ${\mathbf{P}}_{i}$ :

## Eq. 22

${\mathbf{P}}_{i}=\mathbf{M}{\mathbf{P}}_{j},$
or
$\mid \begin{array}{c}{\mathbf{P}}_{1}\\ {\mathbf{P}}_{2}\\ \dots \\ {\mathbf{P}}_{{N}_{i}}\end{array}\mid =\mid \begin{array}{cccccc}{\mathbf{M}}_{1,1}& \dots & {\mathbf{M}}_{{N}_{O},1}& {\mathbf{M}}_{{N}_{O}+1,1}& \dots & {\mathbf{M}}_{{N}_{j},1}\\ {\mathbf{M}}_{1,2}& \dots & {\mathbf{M}}_{{N}_{O},2}& {\mathbf{M}}_{{N}_{O}+1,2}& \dots & {\mathbf{M}}_{{N}_{j},2}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ {\mathbf{M}}_{1,{N}_{i}}& \dots & {\mathbf{M}}_{{N}_{O},{N}_{i}}& {\mathbf{M}}_{{N}_{O}+1,{N}_{i}}& \dots & {\mathbf{M}}_{{N}_{j},{N}_{i}}\end{array}\mid \mid \begin{array}{c}{\mathbf{P}}_{1}\\ \dots \\ {\mathbf{P}}_{{N}_{0}}\\ {\mathbf{P}}_{{N}_{0}+1}\\ \dots \\ {\mathbf{P}}_{{N}_{j}}\end{array}\mid .$
The values of the real sources are constant. Initially, the values of the virtual sources are set to zero. A new set of values for the virtual sources ${\mathbf{P}}_{i}$ is calculated by ${\mathbf{P}}_{i}=\mathbf{M}{\mathbf{P}}_{j}$ . This first calculation of ${\mathbf{P}}_{i}$ is called the Born approximation. The values of ${\mathbf{P}}_{i}$ are used to update the virtual source values in ${\mathbf{P}}_{j}$ :

## Eq. 23

${\mathbf{P}}_{j}\left(j={N}_{0}+1:{N}_{j}\right)={\mathbf{P}}_{i}\left(i=1:{N}_{i}\right).$
Now the process of sequentially using Eqs. 22, 23 is iterated, until the values of ${\mathbf{P}}_{i}$ no longer change. This occurs rapidly for an absorbing object, usually less than nine iterations. For a scattering object (see later), convergence usually takes more iterations, perhaps 20 to 50, but is still rapid. For scattering objects, it is usually very helpful to let each update be a partial update, for example only a $1∕3$ step toward the calculated update value, in other words:

## Eq. 24

${\mathbf{P}}_{j}\left(j={N}_{0}+1:{N}_{j}\right)={\mathbf{P}}_{i}\left(i=1:{N}_{i}\right)+\left[\mathbf{M}{\mathbf{P}}_{j}-i\left(i=1:{N}_{i}\right)\right]∕3.$
This makes convergence smoother and more robust.

If an object is rather large, then it is better to subdivide the object volume into smaller subvolumes, each acting as an independent virtual source. For a scattering object, virtual sources, some positive and some negative, are placed on the surface of the object to create an effective dipole that mimics the behavior of the scattering object. The following discussion outlines these procedures. Figure 8 illustrates the virtual sources described, Fig. 8a showing the absorbing subvolumes and Fig. 8b showing the surface scattering sources that yield a dipole with positive sources on the hemisphere facing upstream into the gradient of fluence rate toward the real source, and negative sources on the hemisphere facing downstream down the gradient of fluence rate away from the real source.

## Fig. 8

The virtual sources for perturbation by an object. (a) Subvolumes within an object, accounting for object absorption and scattering-enhanced absorption. (b) Surface sources on the object’s surface, creating a dipole oriented along the gradient of fluence rate. There are two kinds of perturbation that contribute to the source vector $\mathbf{S}$ and the elements in the $\mathbf{M}$ matrix. One is due to absorption and one is due to scattering. Consider a spherical object of radius $R$ that is both absorbing and scattering. The volume of the object can be subdivided into a set of ${N}_{a}$ subvolumes that fill the object’s volume, and each is called a “virtual absorbing volume source.” The vector ${P}_{j}\left(j={N}_{0}+1:{N}_{0}+{N}_{a}\right)$ is filled with these virtual absorbing volume sources, in this example expressed as small spheres, each of radius $a$ and subvolume ${V}_{i}=\left(4∕3\right)\pi {a}^{3}$ . The radius $a$ is usually chosen to approximately equal ${\mathit{MFP}}^{\prime }∕2$ , which sets the number of small spheres required to fill the object’s volume ${N}_{a}$ , but the total volume of all the small spheres must equal the volume of the object $V={N}_{a}\left(4∕3\right)\pi {a}^{3}$ , so $a={\left[3V∕\left(4\pi {N}_{a}\right)\right]}^{1∕3}$ . The elements of the $\mathbf{M}$ matrix for $i=1:{N}_{a}$ are specified for each $j$ ’th source by:

## Eq. 25

$\mathbf{M}\left(j,i\right)=-\left[\Delta {\mu }_{a}\left(i\right)+\Delta {\mu }_{s}^{\prime }\left(i\right)\frac{{\mu }_{ao}}{{\mu }_{so}^{\prime }}\right]V\left(i\right)T\left(j,i\right).$
The factor $-\Delta {\mu }_{ai}{V}_{i}{T}_{j,i}$ denotes the effect of the $j$ ’th source on the $i$ ’th virtual source, and is the same as the introductory example presented earlier [see Eq. 20]. The second term, $-\Delta {\mu }_{si}^{\prime }\left({\mu }_{a0}∕{\mu }_{s0}^{\prime }\right){V}_{i}{T}_{j,i}$ is a scattering-enhanced absorption due to the extra scattering $\Delta {\mu }_{si}^{\prime }$ that causes light to be trapped within the object, thereby extending its path length within the object.

The surface of the object is studded with point sources called “virtual scattering surface sources,” which mimic the scattering of the object due to its $\Delta {\mu }_{si}^{\prime }$ . The surface is divided into a set of ${N}_{s}$ incremental surfaces, each with a surface area $\Delta {A}_{i}$ . The sum of the incremental surfaces equals the total surface area of the object. The vector ${\mathbf{P}}_{j}\left(j={N}_{0}+{N}_{a}+1:{N}_{j}\right)$ , where ${N}_{j}={N}_{o}+{N}_{a}+{N}_{s}$ , is filled with these virtual scattering surface sources. The elements of the $\mathbf{M}$ matrix for $i={N}_{a}+1:{N}_{i}$ , where ${N}_{i}={N}_{a}+{N}_{s}$ , is specified for each $j$ source by:

## Eq. 26

$\mathbf{M}\left(j,i\right)=\frac{\Delta {\mu }_{a}+\Delta {\mu }_{s}^{\prime }}{{\mu }_{so}^{\prime }}\Delta A\left(i\right)\frac{D}{d}\left[{T}_{1}\left(j,i\right)-{T}_{2}\left(j,i\right)\right],$
where ${T}_{1}\left(j,i\right)$ is the transport from the $j$ ’th source to a position $d∕2$ just outside the object at the $i$ ’th virtual source position, and ${T}_{2}\left(j,i\right)$ is the transport from the $j$ ’th source to a position $d∕2$ just inside the object at the $i$ ’th virtual source position. The distance $d$ is usually chosen to equal ${\mathit{MFP}}^{\prime }$ , such that these two sources inside and outside the object’s surface are sufficiently distant from each other to interact using diffusion theory. Positions 1 and 2 are aligned along the normal to the surface. Consequently, the factor $\left(D∕d\right)\left({T}_{1}-{T}_{2}\right)$ is the flux crossing into the object normal to the surface at the incremental area. The strength of the virtual surface scattering source is proportional to this flux. If flux gradient is driving light into the object, as occurs on the upstream side of the object closest to the real source, then the value of ${\mathbf{P}}_{i}$ is positive. If the flux gradient is driving light out of the object, as occurs on the downstream side of the object away from the real source, then the value of ${\mathbf{P}}_{i}$ is negative. The positive and negative surface sources on the object create a dipole that radiates a dipolar perturbation oriented along the gradient of flux. Light scatters upstream toward the true light source, and a shadow is cast downstream behind the object.

In the matrix $\mathbf{M}$ , a special case occurs when $i$ equals $j$ , and the self-perturbation by a virtual source must be evaluated. This term is here called ${\mathbf{M}}_{j,j}$ . For the scattering surface sources, the value of ${\mathbf{M}}_{j,j}$ is zero, because a point source cannot exert a gradient across itself. For the absorbing volume sources of radius $a$ , the value of ${\mathbf{M}}_{j,j}$ is calculated:

## Eq. 27

${\mathbf{M}}_{j,j}={\int }_{r=0}^{a}\Delta {\mu }_{a}\frac{\mathrm{exp}\left(-r∕\delta \right)}{4\pi Dr}2\pi {r}^{2}dr=\frac{\Delta {\mu }_{a}}{{\mu }_{a0}}\left[1-\left(1+\frac{a}{\delta }\right)\mathrm{exp}\left(-a∕\delta \right)\right].$
After iterating the matrix calculation [Eqs. 22, 23], a stable vector ${\mathbf{P}}_{j}$ is achieved. The fluence rate ${\varphi }_{d}$ at any observation point, for example at a detector, can now be calculated:

## Eq. 28

${\varphi }_{d}=\sum _{j=1}^{{N}_{j}}{\mathbf{P}}_{j}{T}_{j\to d},$
where ${T}_{j\to d}$ is the transport from the $j$ ’th source to the point of observation. All the sources ${\mathbf{P}}_{j}$ , both real and virtual, contribute to the observed fluence rate.

Figure 9 shows an example of a $0.4\text{-}\mathrm{cm}$ -diam. spherical object with adding absorber ( $\Delta {\mu }_{a}=1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_{s}^{\prime }=0$ , ${\mu }_{a0}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu }_{s}^{\prime }=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ). The background light field ${\varphi }_{0}$ is shown in Fig. 9a. The light source is located at $\left(x,y,z\right)=\left(-1,0,0\right)$ . The perturbation itself, ${\varphi }_{p}$ , is shown in Fig. 9b. The object casts a spherical zone of negative perturbation around the object, about $-10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-2}$ effect ( $\mathrm{W}∕{\mathrm{cm}}^{2}$ per W of the source) within the object and $-5\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-2}$ effect in the region surrounding the object. Figure 9c shows the fractional perturbation, ${\varphi }_{p}∕{\varphi }_{0}$ [dimensionless], illustrating that the perturbation is about $-0.30$ within the object and cases a shadow of perturbation about $-0.10$ to $-0.20$ behind the object, downstream from the source. The object’s perturbation can significantly affect the light field downstream from the object, but it cannot overcome the strong light field upstream toward the source.

## Fig. 9

An object with increased absorption ( ${\mu }_{a0}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu }_{s{0}^{\prime }}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_{a}=1{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_{s}^{\prime }=0$ ). (a) The background fluence distribution ${\varphi }_{0}\left(x,z\right)$ at $y=0\phantom{\rule{0.3em}{0ex}}\left[\mathrm{W}∕{\mathrm{cm}}^{2}\right]$ , diffusing from the source. (b) The perturbing fluence ${\varphi }_{p}$ due to the object. (c) The fraction perturbation ${\varphi }_{p}∕{\varphi }_{0}$ . The object asserts a spherical zone of negative perturbation, which casts a perturbing shadow downstream but has little effect upstream toward the source of light. Figure 10 shows the same spherical object, this time with added scattering ( $\Delta {\mu }_{s}^{\prime }=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_{a}=0$ , ${\mu }_{a0}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu }_{s}^{\prime }=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ). This time, the object backscatters a positive perturbation upstream toward the source and casts a negative shadow downstream behind the object. Figure 10c shows that the fractional perturbation ${\varphi }_{p}∕{\varphi }_{0}$ is a little different than the absorption only example [Fig. 9c], but essentially presents the same behavior of casting a shadow downstream behind the object away from the source. These examples illustrate that both increased absorption and increased scattering cast shadows downstream, but have marginal effect upstream. A detector is more sensitive to the presence of a perturbing object if the detector is placed in the shadow cast by the object.

## Fig. 10

An object with increase scattering ( ${\mu }_{a0}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu }_{s{0}^{\prime }}=10{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_{a}=0$ , $\Delta {\mu }_{s}^{\prime }=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ). (a), (b), and (c) are the same as in Fig. 9. The object asserts a dipole pattern of perturbation, in which light is scattered upstream toward the source, and a shadow is cast downstream. However, the perturbation has little effect upstream toward the source of light. In summary, the perturbation method is rather simple and quick, and provides a good estimate of the effect of a perturbing object on the field of light at some distant point of observation, such as a detector. The perturbation method has been used to interpret the gastric endoscopic images of reflectance to determine the size and depth of a subsurface blood vessel, to aid the decision about how to handle a gastric bleeding site,7 and to assess the ability to measure the oxygen saturation of a neonatal brain while still in the birth canal by delivery of light to and collection of escaping optical flux from the maternal abdomen, to aid the decision to deliver the neonate by Caesarean section.8

## 2.5.

### Boundaries

A word on boundaries is worthwhile, since most biophotonic applications involve delivery and/or collection of light through an air/tissue surface. The boundary condition is modeled using the method of images adapted to solve the relationship between fluence rate at the surface and flux escaping the surface. The air/tissue boundary is replaced by an infinite tissue with a source of negative light located above (or outside) the original tissue boundary, which serves to draw light out of the original tissue region, thereby mimicking the loss of light due to escape from the tissue.

Figure 11a schematically depicts light diffusing toward an air/tissue surface, labeled as a $J$ flux, and diffusing away from the surface, labeled as an $I$ flux. The light is shown as a zig-zag trajectory, at $60\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ relative to the $z$ axis, because on average each photon travels a distance $L=2\Delta z$ for every $\Delta z$ movement along the $z$ axis. Consider a test point just below the surface within the tissue, $z=0$ . The fluence rate $\varphi \left(z=0\right)$ at this test point will be twice the value of flux passing along the $z$ axis contributed by both the $I$ and $J$ fluxes, in other words $\varphi \left(0\right)=2\left(I+J\right)$ . As the light hits the surface, there is internal reflectance ${r}_{i}$ of about half the light $\left({r}_{i}\approx 0.50\right)$ . So the value of $I$ is ${r}_{i}J$ . Hence, $\varphi \left(0\right)=2\left(I+J\right)=2\left({r}_{i}J+J\right)=2\left(1+{r}_{i}\right)J$ . The escaping flux is the observable reflectance $R=\left(1-{r}_{i}\right)J$ . Hence, $J=R∕\left(1-{r}_{i}\right)$ . Therefore, the fluence rate at the surface is

## Eq. 29

$\varphi \left(0\right)=2\left(I+J\right)=2R\frac{1+{r}_{i}}{1-{r}_{i}}.$
This is the first key equation. Next, consider Fick’s first law of diffusion. The escaping flux is $R=-D\left(\partial \varphi ∕\partial z\right)$ evaluated at $z\approx 0$ . However, the boundary condition is mimicked by assuming an infinite medium and placing a surface of symmetry outside the true tissue surface at a position $z={z}_{b}$ , where ${z}_{b}$ is a negative value. For each true source of light within the tissue, one places a negative image source at a position symmetrically opposite the true source, outside the original tissue above this plane of symmetry. For example, if a $1\text{-}\mathrm{W}$ point source of light is located at a depth ${z}_{s}$ within the tissue, then a negative point source of equal power is placed at a position ${z}_{i}=2{z}_{b}-{z}_{s}$ .

## Fig. 11

Schematic explanation of mismatched boundary condition. (a) The photons diffuse in a J flux toward the surface and an $I$ flux away from the surface, and on average travel a distance $L=2\Delta z$ per $\Delta z$ movement along the $z$ axis. The internal reflectance ${r}_{i}$ at the surface yields $I={r}_{i}J$ . The escaping reflectance $R=\left(1-{r}_{i}\right)J$ . Hence, the fluence rate at the surface is $\varphi \left(0\right)=2R\left(1+{r}_{i}\right)∕\left(1-{r}_{i}\right)$ , Eq. 29. (b) The 1-D fluence rate distribution due to a planar source at $z={z}_{s}$ and an air/tissue boundary at ${z}_{0}$ . A surface of symmetry is placed as ${z}_{b}$ . If a true positive source of light is present at ${z}_{s}$ within the tissue, a negative image source is placed at ${z}_{i}=2{z}_{b}-{z}_{s}$ outside the tissue. The gradient of light between ${z}_{0}$ and ${z}_{b}$ drives a flux out of the surface, $R=D\varphi \left({z}_{0}\right)∕\mid {z}_{b}\mid$ , Eq. 30. Combining Eqs. 29, 30 yields the position ${z}_{b}=-2\left(1+{r}_{i}\right)∕\left(1-{r}_{i}\right)D$ , which is Eq. 31. The example situation illustrated in Fig. 11b uses 1-D planar diffusion theory ${T}_{\text{planar}}\left(z\right)=\mathrm{exp}\left(-z∕\delta \right)∕\left(2{\mu }_{\alpha }\delta \right)$ to better illustrate the gradient of $\varphi$ between a real source at ${z}_{s}$ and a negative image source at ${z}_{i}$ . When 3-D diffusion theory is used to consider the response to a point source that is close to the surface, the $\varphi$ gradient in the region near the point of light delivery may not be perfectly linear, and the solution may not be accurate in this region. However, the solution yields an accurate solution distant $\left(1\phantom{\rule{0.3em}{0ex}}{\mathit{MFP}}^{\prime }\right)$ from the point source, and the boundary condition discussed here is applicable to 3-D diffusion from a point source at ${z}_{s}$ .

The choice of ${z}_{b}$ is such that the gradient from the tissue surface to this virtual surface of symmetry is $-\partial \varphi ∕\partial z=\varphi \left(0\right)∕{z}_{b}$ . Therefore, the escaping flux is equal to

## Eq. 30

$R=-D\frac{\partial \varphi }{\partial z}=D\frac{\varphi \left(0\right)}{{z}_{b}}.$
This is the second key equation. Combining Eqs. 29, 30 and eliminating $R$ yields

## Eq. 31

${z}_{b}=-2\frac{1+{r}_{i}}{1-{r}_{i}}D.$
Equation 30 describes the boundary condition. The value of internal reflectance ${r}_{i}$ specifies the value ${z}_{b}$ . The value ${r}_{i}$ can be calculated:

## Eq. 32

${r}_{i}={\int }_{\theta =0}^{\pi ∕2}{R}_{\text{Fresnel}}\left(\theta \right)2\pi \phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\theta \right)d\theta ,$
where ${R}_{\text{Fresnel}}$ is the Fresnel reflectance as a function of the angle of incidence $\theta$ relative to the surface normal and the refractive index mismatch at the boundary. The calculation of ${r}_{i}$ is 0.4722 for an air/water interface $\left({n}_{\text{water}}=1.33\right)$ and ${r}_{i}=0.5278$ for an air/tissue interface $\left({n}_{\text{tissue}}=1.4\right)$ . A very good approximation is

## Eq. 33

${r}_{i}=0.6681+0.0636n+0.7099∕n-1.4399∕{n}^{2},$
where $n$ is ${n}_{\text{tissue}}$ , calculated by Egan and Hilgeman9 as cited by Groenhuis Ferwerda, and Ten Bosch.10 It is very good for $1.0 .

Therefore, for a generic air/tissue interface, the value of ${r}_{i}$ is $\sim 0.5$ , and the value of $\left(1+{r}_{i}\right)∕\left(1-{r}_{i}\right)$ is $\sim 3$ . Therefore, ${z}_{b}\approx -\left(2\right)\left(3\right)\left({z}_{s}∕3\right)\approx -2{z}_{s}$ . Hence, the position of the negative image source is placed at ${z}_{i}\approx -5{z}_{s}$ .

Once the positions ${z}_{s}$ and ${z}_{i}$ have been specified, the superposition of the two sources yields the net fluence rate in the tissue:

## Eq. 34

$\varphi \left(r,z\right)=\frac{\mathrm{exp}\left(-{r}_{1}∕\delta \right)}{4\pi D{r}_{1}}-\frac{\mathrm{exp}\left(-{r}_{2}∕\delta \right)}{4\pi D{r}_{2}},$
where
${r}_{1}={\left[{\left(r-{r}_{s}\right)}^{2}+{\left(z-{z}_{s}\right)}^{2}\right]}^{1∕2},$
${r}_{2}={\left[{\left(r-{r}_{s}\right)}^{2}+{\left(z-{z}_{i}\right)}^{2}\right]}^{1∕2}.$
Finally, if one takes the gradient of the fluence rate due to the real and image sources that is normal to the tissue surface, evaluated at $z=0$ , this gradient will drive a flux out of the tissue that predicts the escaping reflectance as a function of radial position $R\left(r\right)$ :

## Eq. 35

$R\left(r\right)={\phantom{\mid }-D\frac{\partial \varphi }{\partial z}\mid }_{z=0}=\frac{1}{4\pi }{z}_{s}\left(\frac{1}{\delta }+\frac{1}{{r}_{1}}\right)\frac{\mathrm{exp}\left(-\mid {r}_{1}\mid ∕\delta \right)}{4\pi D{r}_{1}^{2}}+\left({z}_{s}+2{z}_{b}\right)\left(\frac{1}{\delta }+\frac{1}{{r}_{2}}\right)\frac{\mathrm{exp}\left(-\mid {r}_{2}\mid ∕\delta \right)}{4\pi D{r}_{2}^{2}}.$
This boundary condition works rather well at locations $>1\phantom{\rule{0.3em}{0ex}}{\mathit{MFP}}^{\prime }$ distant from the source at ${z}_{s}$ .11

A common situation is a narrow collimated beam of light illuminating a tissue. The collimated beam launches light into the tissue, and the light acts as if there is a point source at ${z}_{s}=1\phantom{\rule{0.3em}{0ex}}{\mathit{MFP}}^{\prime }$ . Equations 34, 35 yield good approximations to the fluence rate $\varphi \left(z,r\right)$ and escaping flux $R\left(r\right)$ at locations greater than $1\phantom{\rule{0.3em}{0ex}}{\mathrm{MPF}}^{\prime }$ distant from the origin at $\left(r=0\right)$ and greater than $1\phantom{\rule{0.3em}{0ex}}{\mathrm{MPF}}^{\prime }$ from the source at ( $r=0$ , $z={z}_{s}$ ). Monte Carlo simulations are recommended for regions close to the source.

## 3.1.

### Partial Differential Equation Representation of Diffusion

Solution to the diffusion problem can be accomplished through formalized approaches to solving the partial differential equation itself. This can be done numerically or analytically, and the solutions can be derived from the time-domain equation:

## Eq. 36

$\frac{1}{c}\frac{\partial \varphi \left(r,t\right)}{\partial t}-\nabla \cdot D\nabla \varphi \left(\mathbf{r},t\right)+{\mu }_{a}\varphi \left(\mathbf{r},t\right)=S\left(\mathbf{r},t\right),$
where each term on the left is a loss of photons, and the term on the right-hand side describes the source of light. This can also be the frequency domain version of this equation, obtained by Fourier transforming each term of the previous equation:

## Eq. 37

$-\nabla \cdot D\nabla \varphi \left(\mathbf{r},\omega \right)+\left({\mu }_{a}+\frac{i\omega }{c}\right)\varphi \left(\mathbf{r},\omega \right)=S\left(\mathbf{r},\omega \right).$
In the latter equation, if the frequency is zero, then the equation results in the steady-state or continuous wave diffusion equation:

## Eq. 38

$-\nabla \cdot D\nabla \varphi \left(\mathbf{r}\right)+{\mu }_{a}\varphi \left(\mathbf{r}\right)=S\left(\mathbf{r}\right).$
Analytic solutions for the diffusion equations are readily available and can easily be integrated into numerical algorithms to provide Green’s-function-based discrete sum solutions to more complex problems.11, 12, 13 However, the validity and ease of analytic methods tend to break down where the exterior or interior boundaries become complex, thereby requiring more elaborate region integral approaches. Additionally, if the goal of the problem is the combination of clinical or biological imaging tissue volumes where both the boundaries are complex and the interior property values vary, then the optimal choice may transition away from hybrid analytic/numerical solutions toward purely numerical approaches.

Numerical methods to solve the diffusion equation are numerous, including finite difference, finite element and boundary element approaches.14, 15, 16, 17 The major difference between these approaches is in the formalized approach to solving the partial differential equation; however, the choice also leads to major differences in how the solution must be calculated. From a practical standpoint one of the key differences separating the numerical methods is the shape and origin of the mesh of points on which the region to be simulated is discretized. These differences can be seen in the representative meshes displayed in Fig. 12 .

## Fig. 12

Example meshes are shown for (a) finite difference solver in Cartesian coordiates with adaptive mesh refinement in each of the $x,y,z$ directions; (b) finite element solver using triangular elements in two dimensions or tetrahedral elements in three dimensions; and (c) a boundary element solver where the mesh is just composed of nodes at the boundaries outlining homogeneous regions. Often the trickiest part of solving the solution in these domains is how to deal with the boundary between the diffusing regime, and what is outside this domain. Boundary conditions can make solving for analytic solutions more complex, and also the conditions must be validated against either measurements or a Monte Carlo type of simulation.

The choice of numerical approach is often driven by the experience with the method, as each requires an expert to create the code, yet often applications of the codes are better suited to different applications. In Table 4 , a summary of how the mesh and memory used scales with the size of the problem is shown. Additionally, the key advantages and weaknesses of each method are listed briefly.

## Table 4

The three major numerical methods for solving the diffusion equation are listed with associated advantage and disadvantage interms of how they scale as the problem gets larger. Here, N is the number of nodes (unknowns) in the mesh, and B is the bandwidth of the mesh. The memory and size of FEM and FDM assume that the mesh is banded (in 2-D), a banded solver such as LU decomposition is used,14, 15 and a sparse matrix solver is used in 3-D.

IntegratedanalyticmethodArbitrarySmall memory $N$ Fast,conceptually theeasiest approachIntegration over highlyshaped regions andheterogeneties are difficult
FinitedifferencemethodCartesian2-D: ${B}^{*}N$ 3-D: $N$ 2-D: ${B}^{2*}N$ 3-D: $N$ Ease of constructionand discretization,sparse/bandedmatricesDifficult to handleirregular boundaries
Finite elementmethodArbitrary(triangular elementsin 2-D; tetrahedronsin 3-D)2-D: ${B}^{*}N$ 3-D: $N$ 2-D: ${B}^{2*}N$ 3-D: $N$ Arbitrary shapes,flexible,sparse/bandedmatricesVolume meshing in 3-Dcan be difficult forarbitrary shapes
BoundaryelementmethodArbitrary(line segments in2-D; triangularelements in 3-D) ${N}^{2}$ ${N}^{3}$ Arbitrary shapes,nodes only onboundaryFull matrices;assumption of piece-wise constant regions

## 3.2.1.

#### Finite difference method

Finite difference methods (FDM) approaches are routinely solved on a discretized Cartesian mesh, where the volume is arbitrary. However, since the mesh is Cartesian, then the volume must be discretized finely enough to allow accurate representation of the tissue boundaries. Usually this regular Cartesian grid used in discretization is the major limitation of finite difference methods, in that most complex boundaries cannot be discretized finely enough to represent smooth curved boundaries. The finite difference solution is obtained by solving the numerical molecule for a diffusive term, on the nodes given in the mesh.14, 15 The time-domain diffusion equation has a major limitation in this approach, in that the finite difference approach has large time requirements due to the iterative approach to converging on the solution, and so the added need to iterate each time step makes this method intractable for very large discretizations. Nonetheless, for highly complex partial differential equations, and very large problems, the finite difference approach has always been the standard to initiate test solutions, mainly because of its ease of construction and intuitive approach to discretization. The diffusion equation is simply discretized into points $i,j,k$ for the $x,y,z$ coordinates, and the fluence at each location calculated by a straightforward approximation of the derivative by a finite difference, such that it becomes:

## Eq. 39

${D}_{i,j,k}\left(\frac{{\varphi }_{i+1,j,k}-{\varphi }_{i-1,j,k}}{\Delta {L}_{i}^{2}}+\frac{{\varphi }_{i,j+1,k}-{\varphi }_{i,j-1,k}}{\Delta {L}_{j}^{2}}+\frac{{\varphi }_{i,j,k+1}-{\varphi }_{i,j,k-1}}{\Delta {L}_{k}^{2}}\right)+\left[\frac{\left({D}_{i+1,j,k}-{D}_{i,j,k}\right)\left({\varphi }_{i+1,j,k}-{\varphi }_{i,j,k}\right)}{\Delta {L}_{i}^{2}}+\frac{\left({D}_{i,j+1,k}-{D}_{i,j,k}\right)\left({\varphi }_{i,j+1,k}-{\varphi }_{i,j,k}\right)}{\Delta {L}_{j}^{2}}+\frac{\left({D}_{i,j,k+1}-{D}_{i,j,k}\right)\left({\varphi }_{i,j,k+1}-{\varphi }_{i,j,k}\right)}{\Delta {L}_{k}^{2}}\right]+{\mu }_{a}{\varphi }_{i,j,k}={S}_{ijk}.$
This equation assumes a forward differencing approach (note that other differencing methods can be assumed,7 and that the optical properties $D$ and ${\mu }_{a}$ vary spatially and so must be discretized as well as the fluence). This finite difference molecule can then be solved sequentially at each node, and by iterating this over the entire mesh until there is little change between iterations, the solution is reached. This is the most intuitive manner to solve this, called the Jacobi method. Since the Jacobi finite difference method is plagued by slow convergence issues, many alternative methods have been produced to speed up the convergence. Most notable is the implicit solution, where the problem is solved as a matrix inversion, rather than iterative relaxation approach. This requires developing a full matrix representation of how each node in the mesh affects all other nodes and summing up these contributions numerically into a matrix that can be solved. Then the solution is reached simply by inverting this large matrix. This approach obviously has large memory demands for reasonably complex problems, and in contrast the explicit solver requires little memory. Several hybrid approaches have been used, with one notable one being called the Crank-Nicholson alternating difference implicit method. This approach cycles between explicit steps and implicit solution steps to speed up the convergence of the solution, but requires more work to create the initial solution. Most finite difference solvers for regular equations are freely available through academic shareware libraries on the Internet.

Multigrid solver methods have been developed as more elaborate approaches to solving the finite difference method faster.18 The central concept in these approaches is to solve the relaxation on a very coarse grid, and then allow rapid extrapolation of the solution from a coarse grid to a finer grid. The optimally designed approach will actually use several difference mesh grid resolutions, and start out solving the most course grid, then systematically solve the solution on further and further refinements of the mesh. Then, accurate solutions will typically require cycling back and forth on the course grid and fine grid a few times to allow rapid convergence with maximal efficiency. A uniquely high speed solver was developed by the researchers at National Oceanic and Atmospheric Administration (NOAA) and used for diffuse spectroscopy for more than a decade.19, 20, 21

## 3.2.2.

#### Finite element method

The FEM approach to solve partial differential equations is based on solving the equation as integrated on arbitrary discrete elements, called basis functions ${\phi }_{i}$ , which map pieces of the entire domain between discrete points, called nodes,15, 16 as shown in Fig. 12b. This is often called the Galerkin method, and is developed in a formalized approach to discretizing any partial differential equation, resulting in a linear system of equations, which can then be solved by matrix methods. The diffusion equation is usually solved in this approach, in what is called the weak formulation, which is a theoretical framework simplifying the solution found, such that it can be solved in the linear formulation.

In the finite element formalism, the frequency domain solution $\varphi \left(r,\omega \right)$ can be discretized onto the basis functions multiplied by weighting coefficients $\varphi ={\sum }_{j=1}^{N}{\Phi }_{j}{\varphi }_{j}$ , where $\Phi$ is the discrete fluence at each node $j$ with weights $\phi$ , which are linear mapping functions. The discrete values of the fluence are determined as part of the solution process, given simple fixed weight functions. In the Galerkin method of weighted residuals, Eq. 37 is multiplied by an identical set of weighting functions ${\varphi }_{i}$ and integrated over the entire problem domain to give:

## Eq. 40

$⟨-\nabla \cdot D\nabla \Phi {\phi }_{i}⟩+⟨\left({\mu }_{a}+\frac{i\omega }{c}\right)\Phi {\phi }_{i}⟩=⟨S{\phi }_{i}⟩,$
where the notation $⟨uv⟩$ simply denotes integration, such that $⟨uv⟩={\int }_{\Omega }uvdx$ , and $\Omega$ is the entire region of the imaging field over which the integrand is calculated. To avoid the need for second order differentiation of the basis functions, this equation is manipulated further through Green’s theorem to arrive at a form without second derivatives:

## Eq. 41

$⟨D\nabla \Phi \cdot \nabla {\phi }_{i}⟩-\oint \stackrel{̂}{s}\cdot D\nabla \Phi {\phi }_{i}ds+⟨\left({\mu }_{a}+\frac{i\omega }{c}\right)\Phi {\phi }_{i}⟩=⟨S{\phi }_{i}⟩.$
Here $\stackrel{̂}{s}$ is the surface normal from the surface being integrated by $\int$ $ds$ . This representation has two advantages: (1) the second derivative terms have been eliminated, and (2) the surface integral provides a vehicle for implementing the boundary conditions. In fact, its integrand is exactly the normal component of the flux through the boundary surface. At all node points within the finite element mesh that are not on its outer boundary, this term will be zero. In solutions of the prior equation, the integrand in this flux term has been typically approximated by:

## Eq. 42

$\stackrel{̂}{s}\cdot D\nabla \varphi =\alpha \varphi ,$
with the proportionality factor $\alpha$ therefore defined as:

## Eq. 43

$\alpha =\frac{\stackrel{̂}{s}\cdot D\nabla \varphi }{\varphi }=\frac{D}{\varphi }{\phantom{\mid }\frac{\partial \varphi }{\partial r}\mid }_{r=R},$
where $r$ is the normal direction from the surface, along the direction of $\stackrel{̂}{s}$ , at all boundary locations $R$ . In this derivation, $\alpha$ has a theoretical value near $\alpha \approx 0.177$ , for a tissue-air interface.22 While analytical forms of $\alpha$ can be determined, it must be recognized that this coefficient is, at best, a rough estimate of the relative fluence rate gradient proportionality. In reality, diffusion theory may not match the measurements near the source or at the boundaries, due to the irradiance being highly anisotropic. So with the limitations in the accuracy of diffusion theory at the boundary, this parameter may not have any meaningful physical value that would allow a good match to measured data of photon transport in tissue. In practice, this coupling parameter is often chosen to allow a good match between the forward simulations and measured data from tissue phantoms.

Using the weak form of the earlier equation, and discretizing the diffusion coefficient $\left(D=\sum {D}_{i}{\phi }_{i}\right)$ and the absorption coefficient $\left({\mu }_{a}=\sum {\mu }_{ai}{\phi }_{i}\right)$ in the same manner as was done with the fluence, then the equation becomes:

## Eq. 44

$⟨D\Phi {\phi }_{k}\nabla {\phi }_{j}\cdot \nabla {\phi }_{i}⟩+⟨{\mu }_{a}\Phi {\phi }_{l}{\phi }_{i}⟩=⟨S{\phi }_{i}⟩+\alpha \Phi {\phi }_{i}.$
The indices $i$ , $j$ , and $k$ denote the different weighting functions for fluence and tissue properties. The approach is used to create a numerical matrix map $\mathbf{A}$ of the diffusion connectivity between the basis functions, where each element of $\mathbf{A}$ is given by:

## Eq. 45

${a}_{ij}=⟨-\sum {D}_{i}{\phi }_{k}\nabla {\phi }_{j}\nabla {\phi }_{i}-\sum {\mu }_{ai}{\phi }_{i}{\phi }_{l}⟩,$
where ⟨ ⟩ signifies the integration over the entire field. Typically, the basis functions $\phi$ are linear or nearly linear interpolation functions between neighboring nodes (i.e., linearly varying from 1 at the node to 0 at all neighboring nodes), so that all significant interactions only take place between neighboring nodes in this integration. Thus, the resulting matrix is then a continuous map of the diffusion equation, with values represented at the discrete node points, and with only local connectivity between nodes. The basis functions are used to interpolate the values to all points between the nodes. More complex basis function representations exist, but are often not used in simpler solvers.15

The approach to simplifying this equation is too lengthy to explain here, but is well described in Refs. 16, 17, 23. The Galerkin method results in a forward solution that reduces to a matrix equation representation of the diffusion equation:

## Eq. 46

$\left[\mathbf{A}\right]\mathbf{\Phi }=\mathbf{S},$
where $\mathbf{\Phi }$ is the vector of fluence values at all nodes $N$ and $\mathbf{S}$ is the vector of sources at each node, driving the predicted fluence. This fluence is then solved simply by inverting the connectivity matrix $\mathbf{A}$ :

## Eq. 47

$\mathbf{\Phi }={\left[\mathbf{A}\right]}^{-1}\mathbf{S}.$
The diffusion matrix $\mathbf{A}$ is always a banded matrix with many zero elements, and so this can be solved quite efficiently using a specialized banded matrix solver. They key to efficiency here is then using a mesh that has the lowest bandwidth possible. The bandwidth is simply a measure of the node numbering, and indicates the maximum difference between the node numbers that are attached to one another, via being in the same triangular or tetrahedral element. Bandwidth reduction of a generated mesh is often required for efficient solution with a banded matrix solver. However, in 3-D, sparse matrix storage schemes are used for efficiency, and the computational burden increase is usually formidable. Yet several working solutions are freely available and commonly used. Generally though, mesh resolutions above 20,000 to 40,000 nodes can become too large to solve on a single computational engine, and must be solved with parallel computer processor implementation or via commercially optimized solvers.

The boundary conditions are applied by adding to the discretized version of this equation into the $\left[\mathbf{A}\right]$ $\mathbf{\Phi }=\mathbf{S}$ matrix equation.24, 25, 26, 27 This only has additive values at the boundary nodes, and is equal to zero at all internal nodes. Zero fluence boundary conditions are often used only in the simplest of simulations, and are easily implemented by setting the off-diagonal matrix elements to zero along the rows corresponding to boundary nodes, and then allowing the diagonal value to be unity. Extrapolated boundary conditions are implemented by adding additional nodes to the mesh exterior and having measurement surface points within the mesh at the location of the true boundary.

Once the source and boundary conditions are set up in the $\mathbf{A}$ matrix and $\mathbf{S}$ vector, the result from inversion of $\mathbf{A}$ is a discretized map of the optical fluence pattern $\mathbf{\Phi }$ , which is accurate if the domain has been discretized finely enough. Again, this approach has strict limits on the discretization and can run into memory problems for arbitrarily large sized domains.

## 3.2.3.

#### Boundary element method

The boundary element method (BEM) is a less explored area of numerical methods, but is quite useful for problems where the interior properties are known to be homogeneous. This method provides a unique approach where the exterior shape can be complex, and the fluence at the boundaries can then be estimated given knowledge of the homogeneous interior regions. Because of the homogeneity constraints of this approach, it is best implemented when the exterior shapes of tissues and the boundaries of interior organs are well known a priori from some structural imaging modality, and then the light fluences at the boundaries of these tissues can be calculated. This has been shown to be fundamentally useful in electrical impedance problems,28 and has been extended to optical diffusion problems in shape recovery,29 fluorescence prediction30 and spectral region estimation.31 The type of mesh used is simply the exterior boundary mesh of points to discretize the region, as shown in Fig. 12c.

The boundary element solution is actually quite similar to the method of analytic integration, in that its derivation begins by assuming that there is a Green’s function solution for the partial differential equation to be solved, and that this solution could then be used to solve for an arbitrary solution by integration over the boundary shape. Under homogeneity conditions, the diffusion equation is simplified into the modified Helmholtz equation for which the Green’s function can be found to be:

## Eq. 48

${G}_{i}\left(r,{r}_{i}\right)=\frac{\mathrm{exp}\left(\frac{-{k}_{l}\mid r-{r}_{i}\mid }{\sqrt{{D}_{l}}}\right)}{4\pi {D}_{l}\mid r-{r}_{i}\mid },\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{0.3em}{0ex}}3\text{-}\mathrm{D},$
where $\left[{\mu }_{a}\left(r\right)+i\omega ∕c\right]={k}_{l}^{2}$ , and ${k}_{l}$ is constant in subdomain $l$ .

The final version of this expansion for the diffusion equation leads to:31

## Eq. 49

$\left[\mathbf{A}\right]\mathbf{\Phi }-\left[\mathbf{B}\right]\left\{\mathbf{D}\left(\partial \Phi ∕\partial n\right)\right\}=\mathbf{S},$
where the matrices $\mathbf{A}$ and $\mathbf{B}$ in this equation have elements:

## 50.

${a}_{ij}=\left(\Omega ∕4\pi \right){\varphi }_{i}+\int {D}_{i}\left(\partial {G}_{i}∕\partial n\right){\varphi }_{i}ds,$
${b}_{ij}=\int {G}_{i}{\varphi }_{i}ds,$
$\mathbf{S}=\int {s}_{0}{G}_{i}ds,$
where ${s}_{0}$ is the source distribution $ds$ is the discrete integral over the boundary angular space, and $\Omega$ is the solid angle enclosed by the boundary at node $i$ . The complete context of the BEM equation derivation is beyond the scope of this review, but can be found in the papers referenced earlier. The BEM provides a major benefit in obtaining the discretization of the domain, especially in 3-D by requiring only surface meshes. Its main weakness is the use of full matrices, which are required, as compared to the banded or sparse matrices used for FEM calculations.

## 4.1.

### Estimation of Tissue Properties

There are several classes of image or property recovery problems that are solved with these forward diffusion models. The commonality in them is that usually there is some measured dataset ${y}_{d}$ , which is a function of the fluence at the boundary $\varphi$ . This can be simulated by a diffusion theory calculation ${y}_{d}=f\left[\varphi \left({\mu }_{a},{\mu }_{s}^{∕}\right)\right]$ , such as the phase, logarithm of the intensity, relative intensity, intensity ratio, mean transit time, etc. These measurements are then used to predict the properties within some unknown region of the tissue. The terminology commonly used is that the forward problem involves obtaining the fluence given a known set of optical properties. The problem of predicting the optical properties given some measure of the fluence is then called the inverse problem. In the inverse problem, regions to be recovered could be:

• 1. the entire region, assumed to be homogeneous

• 2. layers of tissue

• 3. locally embedded regions within the bulk

• 4. recovery of properties at every point in the tissue.

Much effort has been devoted to all of these problems, and the main issue is one of matrix inversion given the dataset. One common approach to solving this inverse problem is to consider this estimation as a Newton-type problem, where the solution is found by nonlinear iteration to the forward model. The forward model can be analytic or numerical in this case, and the reason for iteration is because the forward diffusion problem is not linear and as such cannot lead to a direct matrix solution for estimating optical properties. Exact linear inversions exist for infinite media situations, and these can be extremely useful as initial guess estimates of bulk tissue.32, 33 However, inversion with models that are not for infinite media do not lead to an explicit inversion formula.

The mathematical formulation would be as follows. Assuming that the initial guess is close to the true vector of property values, the Newton method would assume that the goal is to minimize the square difference between the measured and calculated datasets:

## Eq. 51

$\text{minimize}\phantom{\rule{0.3em}{0ex}}f\left({\mu }_{a},{\mu }_{s}^{∕}\right)={\left\{y-{y}_{d}\left({\mu }_{a},{\mu }_{s}^{∕}\right)\right\}}^{2}.$
Taking the derivative of this and setting it equal to zero leads to:

## Eq. 52

$\mathbf{J}\left\{y-{y}_{d}\left({\mu }_{a},{\mu }_{s}^{∕}\right)\right\}=0,$
where $\mathbf{J}$ is the matrix with elements defined by the derivative of $y$ at each measurement point with respect to the property values ${\mu }_{a}$ or ${\mu }_{s}^{∕}$ at each node, giving matrix elements:

## Eq. 53

${j}_{ik}=\left\{\partial {y}_{di}∕d{\mu }_{ak};\partial {y}_{di}∕d{\mu }_{sk}^{∕}\right\}.$
The way to turn this into an iterative problem is to instead of equating this to zero, rather use a Taylor’s expansion about some initial value set, $\left\{{\mu }_{a0};{\mu }_{s}^{∕}0\right\}$ , such that:

## Eq. 54

$0=\left\{y-{y}_{d}\left({\mu }_{a0},{\mu }_{s0}^{∕}\right)\right\}+\mathbf{J}\left\{\Delta {\mu }_{a},\Delta {\mu }_{s}^{∕}\right\}+\cdots \phantom{\rule{0.2em}{0ex}},$
which leads to the update equation:

## Eq. 55

$\mathbf{J}\left\{\Delta {\mu }_{a},\Delta {\mu }_{s}^{∕}\right\}=\left\{y-{y}_{d}\left({\mu }_{ai},{\mu }_{si}^{∕}\right)\right\}.$
That is solved by inversion of $\left[\mathbf{J}\right]$ ; however, since $\mathbf{J}$ is rarely a square matrix (i.e., the number of nodes is almost never exactly equal to the number of measurements), then the update equation becomes:

## Eq. 56

$\left\{\Delta {\mu }_{a,I},\Delta {\mu }_{si}^{∕}\right\}={\left[{\mathbf{J}}^{T}\mathbf{J}\right]}^{-1}{\mathbf{J}}^{T}\left\{y-{y}_{d}\left({\mu }_{ai},{\mu }_{si}^{∕}\right)\right\}.$
The standard Levenberg-Marquardt rationale for inverting the matrix ${\mathbf{J}}^{T}\mathbf{J}$ is to allow it to be regularized, using a scalar regularization parameter $\lambda$ , which is either empirically or algebraically derived:

## Eq. 57

$\left\{\Delta {\mu }_{a,I},\Delta {\mu }_{si}^{∕}\right\}={\left[{\mathbf{J}}^{T}\mathbf{J}+\lambda \mathbf{I}\right]}^{-1}\left[{\mathbf{J}}^{T}\right]\left\{y-{y}_{d}\left({\mu }_{ai},{\mu }_{si}^{∕}\right)\right\}.$
The theoretical approaches to solving this inversion are quite well developed and vary considerably; however, this description provides the basic Newton approach. The methods for matrix inversion of ${\left[{\mathbf{J}}^{T}\mathbf{J}+\lambda I\right]}^{-1}$ largely depend on the size and how well posed the matrix is. In the smallest case, this could be a single region problem with multiple measurements $NM$ , leading to a $\mathbf{J}$ matrix that is $1×NM$ , and readily solved without regularization. In situations where the regions are in layers or where there are “lumps” to be recovered, then the matrix is still small, and may or may not require regularization. An example of this is shown in Fig. 13 , where the interior of a rat cranium is being imaged with NIR light. Figure 13a shows the MRI image of the cranium, where the brain is visible at the bottom as a lighter gray region, the skull is seen as dark black regions, and the muscle is seen as darker gray through the cranium. Figure 13b is a FEM mesh created to allow diffusion simulation of this tissue, separating the muscle and brain into two distinct regions. The diffuse reconstructed values of the tissue absorption coefficient are shown in Fig. 13d where the differences in absorption between brain and muscle tissue are apparent, even in this ill-posed diffuse tomography image. In the full image reconstruction problem, where there are $NN$ nodes and $NM$ measurements, the matrix $\mathbf{J}$ is $NN×NM$ and can require considerable effort to regularize and invert, and still will lead to diffuse blurry images, of the type shown in Fig. 13d. Diffuse tomography always leads to these low resolution types of images, because the forward diffusion transport process makes the image reconstruction problem highly ill-posed and nonunique. Virtually all diffuse tomography imaging is limited in this manner, in that the abilty to recover sharp features is significantly diminished. Resent work in NIR imaging has been in the area of using diffuse spectroscopy as an adjuvant to standard imaging systems, such as MRI or CT. An example of how this would be done is shown in Fig. 13c, where the same rodent brain as before is recovered as a predefined region in the absorption coefficient. In this approach, it is possible to obtain quantitatively accurate bulk tissue values, because the inversion problem is no longer one dominated by image reconstruction, since the recovery only involves two regions instead of $NN$ regions. It has been shown that image-guided recovery of NIR measurements allows much more accurate spectroscopy of tissue, if there is a meaningful way to combine NIR with an imaging system to give the prior information.

## Fig. 13

Recovery of absorption coefficient at $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ is shown for a rat brain region,58 as shown by (a) the coronal MRI image. (b) The brain and muscle tissues can be segmented into two material properties, and then the property recovery problem can either be done as only two regions (i.e., Jacobian matrix size $2×NM$ ) as shown in (c), or as a full diffuse image reconstruction (i.e., Jacobian matrix size $NN×NM$ ) as shown in (d). Recovery of tissue values at a single wavelength are readily achieved, and computer codes such as TOAST (University College London34) and NIRFAST (Dartmouth College and University of Exeter35) are freely available for this problem.

## 4.2.

### Absorption Spectroscopy

Diffuse spectroscopy of tissue is often the primary goal in diffuse optics, and the nature of the problem can be transmission spectroscopy for absorption or scatter, or luminescence spectroscopy for bioluminescence, fluorescence, or phosphorescence. Any problem that requires multiple wavelength fitting to a predetermined extinction or absorption spectrum can fall into the secondary inversion problem of spectroscopy.

In absorption spectroscopy, the goal is to then fit the measured or computed absorption coefficients (e.g., from solving the inverse problem in previous section) to the known molar absorption spectra of the features that absorb in tissue, such as hemoglobins, breakdown products of heme, water, lipids, cytochromes, and potentially other absorbers that are endogenous or exogenous. The problem is formulated as:

## Eq. 58

${\mu }_{a}=\left[\mathbf{\epsilon }\right]\mathbf{C},$
where spectral fitting is ideally achieved by simple inversion of matrix $\epsilon$ , with elements of molar absorption for each wavelength. The inversion is then:

## Eq. 59

$\mathbf{C}={\left[{\mathbf{\epsilon }}^{T}\mathbf{\epsilon }\right]}^{-1}{\epsilon }^{T}{\mu }_{a}.$
Since the absorption relationship to concentration is linear, this is a well-posed problem with little need for regularization. In the regular tissue spectroscopy problem, the number of absorbers may be less than six, and the number of measurements can be high, leading to an overdetermined problem, although noise in the data is known to lead to inaccurate measurements. Most of the major endogenous absorbers of soft tissue in the diffuse transport regime are shown in Fig. 14a , to show the trough in absorption between hemoglobin and water, bounding the diffuse regime on the red and NIR regions, respectively. Outside the range of $600\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1100\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , it is difficult to believe that diffusion theory would be an acceptable approximation, because in many tissues the absorption coefficient is of similar value to the transport scattering coefficient, as seen in this graph.

## Fig. 14

(a) The absorption spectra of the major chromophores in tissue are shown along with a typical reduced scattering spectra. The changes to the reduced scattering spectra with Mie-like scatterer size change are shown in (b).59 The luminescence spectra for several known endogenous emitters is shown in (c). Spectroscopy of tissue with this approach to fitting the spectra has been shown to be enormously successful in breast tumor spectroscopy,36 brain tissue monitoring,37 and muscle physiology,38, 39 to name a few applications. Image-guided diffuse optical spectroscopy has recently been established as a useful tool to image cancer tumors, and uses the principles of region-based estimation to reduce the size of the Jacobian and allow more accurate recovery of the interior regions. Additionally, incorporating the spectral inversion into the image recovery problem has been demonstrated to allow more accurate estimation of the tissue values.40, 41, 42 It is generally assumed that addition of any constraints into the inversion problem, in the way of spectral or spatial prior knowledge, will allow a well-posed inverse problem and potentially more accurate estimation of tissue properties.

## 4.3.

### Elastic Scatter Spectroscopy

In scatter spectroscopy, the nature of the scatter is elastic scattering, and known to fit to Mie theory under certain circumstances. Several research groups have shown the ability to fit the scattered spectrum to the effective particle size and number density, under the assumption of a Mie-like scattering model. Alternatively, the spectra is thought to fit to a Mie-like spectra with some Rayleigh scatter influence, which has been shown to be near:

## Eq. 60

${\mu }_{s}^{∕}\left(\lambda \right)=a{\lambda }^{-b}+c{\lambda }^{-4}.$
Since this is a nonlinear equation, it is again not easily inverted for $a$ or $c$ , and ultimately these parameters are not intuitive, as the values for $a$ and $b$ are arbitrary from an analytic approximation to the Mie theory. The important spectra to fit comes from the Mie theory:

## Eq. 61

${\mu }_{s}^{\prime }\left(\lambda \right)={N}_{o}\sum _{i=1}^{p}f\left({a}_{i}\right)\left(\pi {a}_{i}^{2}\right){Q}_{\mathrm{scat}}\left(m,{a}_{i},\lambda \right)\left[1-g\left(m,{a}_{i},\lambda \right)\right],$
where $f\left({a}_{i}\right)$ is the particle size histogram, with diameter ${a}_{i}$ and index ratio $m$ , with scattering efficiency ${Q}_{\mathit{scat}}$ from Mie theory43 and calculated anisotropy parameter $g$ .44 An illustration of the change in the scattering spectrum with size of the scatterer is shown in Fig. 14b.

Fitting remitted or transmitted scatter spectra to models of the scatterer particle size has become an extremely active area of research, with extensive preclinical and clinical data indicating that scattering spectra can be used to estimate scatterer composition and size in certain cases.45, 46, 47, 48, 49, 50

## 4.4.

### Luminescence Spectroscopy: Fluorescence, Phosphorescence, and Bioluminescence

Fitting of luminescence spectroscopy is less common, but becoming increasingly important for recovery of a unique concentration value,51 and for applications where the spectral shape is altered on remission from the tissue, providing some measure of depth from which the signal is coming.52 This fitting is identical to the absorption spectra fitting, although often there is need for addition of a background spectrum or contaminant spectra that must be included to allow accurate fitting. Examples of fluorescence spectra from relevant fluorophores, which are used routinely in tissue imaging, are shown in Fig. 14c. The amplitudes of each are normalized to unity at their maximum, but the yield of fluorescence is directly proportional to (1) the molar extinction coefficient at the excitation light wavelength, (2) the intensity of the excitation light, (3) the fluorescence quantum yield of the fluorophore within tissue, and (4) the concentration of the fluorophore present. A reduction in any of these four things will lead to reduced fluorescence from the tissue sampled.

Fluorescence tomography for small animal imaging has been largely established by Zacharakis and by Ntziachristos 53, 54, 55, 56 in recent years. Yet the majority of small animal luminescence imaging with commercial instrumentation is done without much consideration of diffusion issues or compensation for diffusion, which can suffice for relative comparisons of fluorophore production or uptake. However, when the location of the fluorophore changes, or when the tissue shape change is relevant, then more extensive modeling is required. Additionally, it has become well established in recent years that if surface measurements are sufficient, then the reflectance geometry of having the source and detector on the same side of the tissue is good enough for imaging. However, if the fluorophore is deeper in the tissue, then the transmission geometry provides a superior signal-to-background value.57 This transmission-based measurement for uptake of fluorophores in tumor tissue will become more utilized as this mode of measurement becomes standard in commercial small animal imaging systems.

## Conclusion

This work focuses on the fundamental concepts leading to the mathematical and numerical estimations of light diffusion in tissue. While the theory can be developed from a number of different perspectives, the focus here is on phenomenological description, allowing physical insights into the problem. Analytic derivations are kept to a minimum while outlining the overall concepts. In the numerical methods presented, the focus is on providing a practical outlines of the strengths and weaknesses of the different methods. This review is a brief introduction to what has become a well developed mature field of diffuse light transport modeling. The references listed are a good introduction to the fundamentals and applications that have emerged in the past few decades of work.

## References

1.

A. J. Welch and M. J. C. van Gemert, Optical-Thermal Response of Laser-Irradiated Tissue, Plenum Press, New York (1995). Google Scholar

2.

L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML-Monte-Carlo modeling of light transport in multilayered tissues,” Comput. Methods Programs Biomed., 47 (2), 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar

3.

L. H. Wang and S. L. Jacques, “Hybrid model of Monte-Carlo simulation and diffusion-theory for light reflectance by turbid media,” J. Opt. Soc. Am. A Opt. Image Sci. Vis, 10 (8), 1746 –1752 (1993). 1084-7529 Google Scholar

5.

L. O. Svaasand, B. J. Tromberg, R. C. Haskell, T. T. Tsay, and M. W. Berns, “Tissue characterization and imaging using photon density waves,” Opt. Eng., 32 (2), 258 –266 (1993). https://doi.org/10.1117/1.2170990 0091-3286 Google Scholar

6.

M. R. Ostermeyer and S. L. Jacques, “Perturbation theory for diffuse light transport in complex biological tissues,” J. Opt. Soc. Am. A Opt. Image Sci. Vis, 14 (1), 255 –261 (1997). 1084-7529 Google Scholar

7.

C. J. Gostout and S. L. Jacques, “Infrared video imaging of subsurface vessels—a feasibility study for the endoscopic management of gastrointestinal-bleeding,” Gastrointest. Endosc., 41 (3), 218 –224 (1995). 0016-5107 Google Scholar

8.

S. L. Jacques, N. Ramanujam, G. Vishnoi, R. Choe, and B. Chance, “Modeling photon transport in transabdominal fetal oximetry,” J. Biomed. Opt., 5 (3), 277 –282 (2000). https://doi.org/10.1117/1.429996 1083-3668 Google Scholar

9.

W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials, Academic, New York (1979). Google Scholar

10.

R. A. J. Groehnuuis, H. A. Ferwerda, and J. J. Ten Bosch, “Scattering and absorption of turbid materials determined from reflectaion measurements. 1: Theory,” Appl. Opt., 22 (16), 2456 –2462 (1983). 0003-6935 Google Scholar

11.

M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical-properties,” Appl. Opt., 28 (12), 2331 –2336 (1989). 0003-6935 Google Scholar

12.

S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol., 37 (7), 1531 –1560 (1992). https://doi.org/10.1088/0031-9155/37/7/005 0031-9155 Google Scholar

13.

B. W. Pogue and M. S. Patterson, “Frequency domain optical absorption spectroscopy of finite tissue volumes using diffusion theory,” Phys. Med. Biol., 39 1157 –1180 (1994). https://doi.org/10.1088/0031-9155/39/7/008 0031-9155 Google Scholar

14.

L. L. Lapidus and G. F. Pinder, Numerical Solution of Partial Differential Equations in Science and Engineering, John Wiley and Sons, New York (1999). Google Scholar

15.

D. R. Lynch, Numerical Partial Differential Equations for Environmental Scientists and Engineers—A First Practical Course, Springer, New York (2005). Google Scholar

16.

K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys., 22 (6), 691 –701 (1995). https://doi.org/10.1118/1.597488 0094-2405 Google Scholar

17.

S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite-element approach for modeling photon transport in tissue,” Med. Phys., 20 (2), 299 –309 (1993). https://doi.org/10.1118/1.597069 0094-2405 Google Scholar

18.

B. W. Pogue, M. S. Patterson, H. Jiang, and K. D. Paulsen, “Initial assessment of a simple system for frequency-domain diffuse optical tomography,” Phys. Med. Biol., 40 (10), 1709 –1729 (1995). https://doi.org/10.1088/0031-9155/40/10/011 0031-9155 Google Scholar

19.

J. C. Adams, “MUDPACK: multigrid portable fortran software for the efficient solution of linear elliptic differential equations,” J. Appl. Math. Comput., 34 I13 –46 (1989). Google Scholar

20.

B. W. Pogue, M. S. Patterson, H. Jiang, and K. D. Paulsen, “Initial assessment of a simple system for frequency domain diffuse optical tomography,” Phys. Med. Biol., 40 1709 –1729 (1995). https://doi.org/10.1088/0031-9155/40/10/011 0031-9155 Google Scholar

21.

D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick–Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media,” Appl. Opt., 36 2260 –2272 (1997). 0003-6935 Google Scholar

22.

B. W. Pogue, S. Geimer, T. O. McBride, S. Jiang, U. L. Österberg, and K. D. Paulsen, “Three-dimensional simulation of near-infrared diffusion in tissue: boundary condition and geometry analysis for finite element image reconstruction,” Appl. Opt., 40 (4), 588 –600 (2001). https://doi.org/10.1364/AO.40.000588 0003-6935 Google Scholar

23.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 (2), R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar

24.

R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A, 11 (10), 2727 –2741 (1994). 0740-3232 Google Scholar

25.

R. Aronson, “Boundary-conditions for diffusion of light,” J. Opt. Soc. Am. A, 12 (11), 2532 –2539 (1995). 0740-3232 Google Scholar

26.

A. H. Hielscher, L. Wang, K. Tittel, and S. Jacques, “Influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol., 40 1957 –1975 (1995). https://doi.org/10.1088/0031-9155/40/11/013 0031-9155 Google Scholar

27.

M. S. Patterson and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt., 28 2331 –2336 (1989). 0003-6935 Google Scholar

28.

S. Babaeizadeh, D. H. Brooks, D. Isaacson, and J. C. Newell, “Electrode boundary conditions and experimental validation for BEM-based EIT forward and inverse solutions,” IEEE Trans. Med. Imaging, 25 (9), 1180 –1188 (2006). https://doi.org/10.1109/TMI.2006.879957 0278-0062 Google Scholar

29.

J. Sikora, A. Zacharopoulos, A. Douiri, M. Schweiger, L. A. Horesh, and J. Ripoll, “Diffuse photon propagation in multilayered geometries,” Phys. Med. Biol., 51 497 –516 (2006). https://doi.org/10.1088/0031-9155/51/3/003 0031-9155 Google Scholar

30.

F. Fedele, J. P. Laiblea, and M. J. Eppstein, “Fluorescence photon migration by the boundary element method,” J. Comp. Physiol., 210 109 –132 (2005). 0373-0859 Google Scholar

31.

S. Srinivasan, B. W. Pogue, C. Carpenter, P. Yalavarthy, and K. D. Paulsen, “A boundary element approach for image-guided near-infrared absorption and scatter estimation,” Med. Phys., 0094-2405 Google Scholar

32.

S. Fantini, M. A. Franceschini, and E. Gratton, “Semi-infinite-geometry boundary-problem for light migration in highly scattering media—a frequency-domain study in the diffusion-approximation,” J. Opt. Soc. Am. B, 11 (10), 2128 –2138 (1994). 0740-3224 Google Scholar

33.

B. W. Pogue, K. D. Paulsen, H. Kaufman, and C. Abele, “Calibration of near infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating phantoms,” J. Biomed. Opt., 5 (2), 182 –93 (2000). 1083-3668 Google Scholar

36.

B. J. Tromberg, A. Cerussi, N. Shah, M. Compton, A. Durkin, D. Hsiang, J. Butler, and R. Mehta, “Imaging in breast cancer: diffuse optics in breast cancer: detecting tumors in pre-menopausal women and monitoring neoadjuvant chemotherapy,” Breast Cancer Res. Treat., 7 (6), 279 –285 (2005). 0167-6806 Google Scholar

37.

D. A. Boas, A. M. Dale, and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” Neuroimage, 23 S275 –S288 (2004). 1053-8119 Google Scholar

38.

E. Gratton, S. Fantini, M. A. Franceschini, G. Gratton, and M. Fabiani, “Measurements of scattering and absorption changes in muscle and brain,” Philos. Trans. R. Soc. London, Ser. B, 352 (1354), 727 –735 (1997). 0962-8436 Google Scholar

39.

T. Hamaoka, T. Katsumura, N. Murase, S. Nishio, T. Osada, T. Sako, H. Higuchi, Y. Kurosawa, T. Shimomitsu, M. Miwa, and B. Chance, “Quantification of ischemic muscle deoxygenation by near infrared time-resolved spectroscopy,” J. Biomed. Opt., 5 (1), 102 –105 (2000). https://doi.org/10.1117/1.429975 1083-3668 Google Scholar

40.

A. Li, G. Boverman, Y. Zhang, D. Brooks, E. L. Miller, M. E. Kilmer, Q. Zhang, E. M. C. Hillman, and D. A. Boas, “Optimal linear inverse solution with multiple priors in diffuse optical tomography,” Appl. Opt., 44 (10), 1948 –1956 (2005). https://doi.org/10.1364/AO.44.001948 0003-6935 Google Scholar

41.

S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, J. J. Gibson, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “In vivo hemoglobin and water concentration, oxygen saturation, and scattering estimates from near-infrared breast tomography using spectral reconstruction,” Acad. Radiol., 1076-6332 Google Scholar

42.

B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, T. Tosteson, J. B. Weaver, S. P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures using hybrid MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. U.S.A., 103 (23), 8828 –8833 (2006). https://doi.org/10.1073/pnas.0509636103 0027-8424 Google Scholar

43.

C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley and Sons, New York (1983). Google Scholar

44.

H. J. van Staveren, C. J. M. Moes, J. van Marle, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in intralipid-10% in the wavelength range of $400–1100\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$,” Appl. Opt., 30 (31), 4507 –4514 (1991). 0003-6935 Google Scholar

45.

V. Backman, M. B. Wallace, L. T. Perelman, J. T. Arendt, R. Gurjar, M. G. Muller, Q. Zhang, G. Zonios, E. Kline, J. A. McGilligan, S. Shapshay, T. Valdez, K. Badizadegan, J. M. Crawford, M. Fitzmaurice, S. Kabani, H. S. Levin, M. Seiler, R. R. Dasari, I. Itzkan, J. Van Dam, M. S. Feld, and T. McGillican, “Detection of preinvasive cancer cells,” Nature (London), 406 (6791), 35 –36 (2000). https://doi.org/10.1038/35017638 0028-0836 Google Scholar

46.

R. S. Gurjar, V. Backman, L. T. Perelman, I. Georgakoudi, K. Badizadegan, I. Itzkan, R. R. Dasari, and M. S. Feld, “Imaging human epithelial properties with polarized light-scattering spectroscopy,” Nat. Med., 7 (11), 1245 –1248 (2001). 1078-8956 Google Scholar

47.

A. Wax, C. H. Yang, M. G. Muller, R. Nines, C. W. Boone, V. E. Steele, G. D. Stoner, R. R. Dasari, and M. S. Feld, “In situ detection of neoplastic transformation and chemopreventive effects in rat esophagus epithelium using angle-resolved low-coherence interferometry,” Cancer Res., 63 (13), 3556 –3559 (2003). 0008-5472 Google Scholar

48.

A. Dhar, K. Johnson, M. Novelli, S. Bown, I. Bigio, L. Lovat, and S. Bloom, “Elastic scattering spectroscopy for the diagnosis of colonic lesions: initial results of a novel optical biopsy technique,” Gastrointest. Endosc., 63 (2), 257 –261 (2006). 0016-5107 Google Scholar

49.

X. Wang, B. W. Pogue, S. Jiang, H. Dehghani, X. Song, S. Srinivasan, B. A. Brooksby, K. D. Paulsen, C. Kogel, A. P. Poplack, and W. A. Wells, “Image reconstruction of effective mie scattering parameters of breast tissue in vivo with near-infrared tomography,” J. Biomed. Opt., 11 (4), 041106 (2006). https://doi.org/10.1117/1.2342747 1083-3668 Google Scholar

50.

C. M. Carpenter, B. W. Pogue, S. Jiang, H. Dehghani, X. Wang, K. D. Paulsen, W. A. Wells, J. Forero, C. Kogel, J. B. Weaver, S. P. Poplack, and P. A. Kaufman, “Image-guided optical spectroscopy provides molecular-specific information in vivo: MRI-guided spectroscopy of breast cancer hemoglobin, water, and scatterer size,” Opt. Lett., 32 (8), 933 –935 (2007). https://doi.org/10.1364/OL.32.000933 0146-9592 Google Scholar

51.

R. M. Levenson and J. R. Mansfield, “Spectral imaging in biology and medicine: slices of life,” Cytom. A., 69 (8), 748 –758 (2006). Google Scholar

52.

H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally-resolved bioluminescence optical tomography,” Opt. Lett., 31 (3), 365 –367 (2006). https://doi.org/10.1364/OL.31.000365 0146-9592 Google Scholar

53.

G. Zacharakis, J. Ripoll, R. Weissleder, and V. Ntziachristos, “Fluorescent protein tomography scanner for small animal imaging,” IEEE Trans. Med. Imaging, 24 (7), 878 –885 (2005). https://doi.org/10.1109/TMI.2004.843254 0278-0062 Google Scholar

54.

V. Ntziachristos, E. A. Schellenberger, J. Ripoll, D. Yessayan, E. Graves, A. Bogdanov Jr., L. Josephson, and R. Weissleder, “Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin V-Cy5.5 conjugate,” Proc. Natl. Acad. Sci. U.S.A., 101 (33), 12294 –12299 (2004). https://doi.org/10.1073/pnas.0401137101 0027-8424 Google Scholar

55.

V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med., 8 (7), 757 –760 (2002). 1078-8956 Google Scholar

56.

V. Ntziachristos, C. Bremer, E. E. Graves, J. Ripoll, and R. Weissleder, “In vivo tomographic imaging of near-infrared fluorescent probes,” Mol. Imaging, 1 (2), 82 –88 (2002). https://doi.org/10.1162/153535002320162732 1535-3508 Google Scholar

57.

V. Ntziachristos, G. Turner, A. Dunham, S. Windsor, A. Soubret, J. Ripoll, and H. A. Shih, “Planar fluorescence imaging using normalized data,” J. Biomed. Opt., 10 064007 (2005). https://doi.org/10.1117/1.2136148 1083-3668 Google Scholar

58.

H. Xu, “MRI-coupled broadband near-infrared tomography for small animal brain studies,” Thayer School of Engineering, 195 (2005) Google Scholar

59.

X. Wang, B. W. Pogue, S. Jiang, X. Song, K. D. Paulsen, C. Kogel, S. P. Poplack, and W. A. Wells, “Approximation of Mie scattering parameters from near-infrared tomography of health breast tissue in vivo,” J. Biomed. Opt., 10 (5), 051704-1--10 (2005). https://doi.org/10.1117/1.2098607 1083-3668 Google Scholar
©(2008) Society of Photo-Optical Instrumentation Engineers (SPIE)
Steven L. Jacques and Brian W. Pogue "Tutorial on diffuse light transport," Journal of Biomedical Optics 13(4), 041302 (1 July 2008). https://doi.org/10.1117/1.2967535
Published: 1 July 2008
JOURNAL ARTICLE
19 PAGES SHARE