Abstract

We reconstruct images of the absorption and the scattering coefficients for diffuse optical tomography using five different models for light propagation in tissues: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation. The last four models listed are approximations of the radiative transport equation that take into account forward-peaked scattering analytically. Using simulated data from the numerical solution of radiative transport equation, we solve the inverse problem for the absorption and scattering coefficients using the transport-backtransport method. Through comparison of the numerical results, we show that all of these light scattering models produce good image reconstructions. In addition, we show that these approximations afford considerable computational savings over solving the radiative transport equation. However, all of the models exhibit significant “crosstalk” between absorption and scattering coefficient images. Among the approximations, we have found that the generalized Fokker-Planck-Eddington equation produced the best image reconstructions in comparison with the image reconstructions produced by the radiative transport equation.

©2009 Optical Society of America

1. Introduction

In diffuse optical tomography, one seeks to reconstructs the optical properties of tissues using measurements of scattered near-infrared light. This imaging modality has applications in neonatal brain imaging and breast imaging, for example (see [1, 2] and references contained therein). The major challenge in this field is taking into account the strong multiple scattering of light in tissues.

Light propagation in tissues is governed by the theory of radiative transport [3]. This theory takes into account absorption and scattering due to inhomogeneities. The radiative transport equation is an integral-differential equation for the specific intensity. This equation is valid when distances between scatterers are much larger than the wavelength. Consequently, this equation does not take into account interference effects between scatterers. There are very few specific cases for which we know the analytical solution to the radiative transport equation. In addition, it is difficult to solve the radiative transport equation numerically, especially when the medium is optically thick or scatters with a sharp peak in the forward direction. Both of these situations are relevant to light propagation in tissues. Forward-peaked scattering requires restrictively large computational resources to approximate the integral operator in the transport equation accurately.

There are several approximations to the radiative transport equation that take into account forward-peaked scattering [4]. Comparisons of the reflectance and transmittance of light by tissues computed using these approximations show good agreement with those computed by the radiative transport equation. These approximations require less computational resources to solve than the radiative transport equation. Hence, these approximations will be useful for diffuse optical tomography if they produce good quality image reconstructions.

To investigate the quality of image reconstructions that these approximations produce, we solve inverse problems using these approximations as models for the measured data and examine the quality of the images they produce. In particular, we compute image reconstructions produced by five light scattering models: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation. Our simulated data is computed using the radiative transport equation with forward-peaked scattering governed by the Henyey-Greenstein scattering phase function.

For the solution of the inverse problem, we use the transport-backtransport method [5, 6]. The transport-backtransport method is one example of the more general propagation-backpropagation method. The propagation back-propagation method has been used widely in medical imaging [1, 7, 8], geophysical applications [9, 10, 11, 12] and other general applications as shown in [13, 14]. This inversion method is comprised of two parts: (1) solving the direct problem and (2) calculating updates for the parameters that we seek to reconstruct. We use the five light scattering models mentioned above for part (1) of the inversion method to reconstruct five different sets of images. Through comparison of these five sets of reconstructed images, we determine which of the approximations are suitable for diffuse optical tomography.

Because we are using the transport-backtransport method to solve the inverse problem, we solve only one direct and one adjoint problem per iteration regardless of the number of unknowns in the problem. In contrast, direct gradient calculation methods require the solution of one direct problem for each parameter in the problem. For the problems we study here, we are dealing with approximately 2500 parameters. Hence, direct gradient calculation methods would require restrictively large computations.

The inverse problem involves the reconstruction of the spatial distribution of the absorption and scattering coefficients from scattered light measurements at the boundary. This problem is severely ill-posed due to multiple scattering. As a result, we do not anticipate a perfect, detailed reconstruction. In fact, one artifact that we will see in our reconstructions is the so-called “crosstalk” in which images of the absorption coefficient show artifacts from the scattering coefficient, and vice versa. As a result, one method to compare the performance of the approximations in image reconstructions is to evaluate the amount of cross-talk in the images.

The remainder of this paper is as follows. In Section 2 we discuss the direct and inverse problems for diffuse optical tomography. Included in that discussion is a description of the radiative transport equation, which governs light propagation in tissues and the transport-backtransport method which we use to solve the inverse problem. In Section 3 We discuss the different approximations of the radiative transport equation for forward-peaked scattering. We give the numerical method we use to solve the direct problem in Section 4. In Section 5 we show the results of the numerical simulations. Section 6 gives our conclusions. Appendix A gives a proof of a result used in the transport-backtransport method in Section 2. Appendix B discusses the self-adjoint properties of the light scattering models we consider here.

2. Diffuse optical tomography

For diffuse optical tomography, we study two mathematical problems: the direct problem and the inverse problem. For the direct problem, we seek to model measurements of scattered light for a given set of optical properties. For the inverse problem, we compare our model for the measurements with the actual measurements to reconstruct the optical properties. In what follows, we discuss both the direct and inverse problems.

2.1. The direct problem

To model light propagation in tissues, we solve the radiative transport equation:

c1tI+ŝ·I+μaIμsI=QonSn1×Ω×[0,T].

Here, S n-1 denotes the set of unit direction vectors, Ω ⊂ Rn denotes the spatial domain with boundary ∂Ω and [0,T] is the interval of time we consider. The absorption coefficient is denoted by μa, the scattering coefficient is denoted by μs and the source is denoted by Q. The scattering operator 𝓛 is defined as

I=I+Sn1f(ŝ·ŝ)I(ŝ,r,t)dŝ.

The scattering phase function f gives the fraction of light scattered in direction ŝ due to light incident in direction ŝ′. We discuss specific choices for 𝓛 in Section 3.

In addition to Eq. (1), we must prescribe boundary and initial conditions. Boundary conditions give the light for all directions pointing into the spatial domain for all points at the boundary over all time. For our boundary conditions, we prescribe that

I=0onΓin={(ŝ,rb,t)Sn1×Ω×[0,T],ν(rb)·ŝ<0},

with ν(r b) denoting the unit outward normal at r b ∈ ∂D. The initial condition gives the intensity at t = 0 for all directions and all points in the spatial domain. For the initial condition, we prescribe that

It=0=0onSn1×Ω.

Hence, the only source of light in the direct problem is given by the source Q. This source Q may correspond to a source at the boundary of the spatial domain.

2.2. The inverse problem

For measurements, we consider angle-averaged quantities at the boundary that take the form

Mμa,μs(rb,t)=Soutn1ŝ.ν(rd)I(ŝ,rd,t)dŝonΩ×[0,T].

Here, S out n-1 = {ŝ ∈ S n-1,ŝ · ν(r d) < 0} with ν(r d) denoting the unit outward normal at r d ∈ ∂Ω, is the set of all angles pointing out of the spatial domain. Suppose we have collected a set of measurements of the form given by Eq. (5) for each source within a set of sources. For the inverse problem, we seek to reconstruct μa and μs using this data and the model for measurements Eq. (5) with I denoting the solution of the direct problem given by Eq. (1) subject to Eq. (3) and Eq. (4).

This inverse problem is nonlinear since the solution of the direct problem depends nonlinearly on the parameters μa and μs. Moreover, this inverse problem is severely ill-posed due to multiple scattering. Finally, we will study inverse problems that are underdetermined since it is typical in practical situations to have only a few detectors placed around the spatial domain.

To address these challenges, we seek to solve the inverse problem by solving an associated optimization problem. We denote the space of μa or μs by P and the space of data M μa,μs by D. The residual is given by the difference between the real data F generated by experimental measurements and the modeled data (i.e. the whole set of sources and receivers) corresponding to the parameter distribution (μa,μs):

R:P×PD,R(μa,μs)=Mμa,μsF.

We introduce the cost functional:

K(μa,μs)=12R(μa,μs)L22.

Now, we turn our attention toward determining the parameter distributions for μa and μs that minimize K. That minimum norm solution yields the image reconstructions for μa and μs that we are seeking.

2.3. The transport-backtransport method

We use the transport-backtransport method developed by O. Dorn [5, 6] to determine the parameter distributions for μa and μs that minimize K given by Eq. (7). We do not describe all of the details of this method here. Rather, we give only the main results.

By perturbing the parameters about μ̅a and μ̅s as μa = μ̅a+δμa and μs = μ̅s+δμs, we find that the linearized residual operator about (μ̅a,μ̅s) is given by

R(μa,μs)R(μa̅,μ̅s)+R'(μa̅,μ̅s)(δμa,δμs),

with R′ denoting the Fréchet derivative of R [15]. By setting the right-hand side of Eq. (8) equal to zero, we obtain a linear system for the perturbations (δμa, δμs). The minimum norm solution of this linear system is given by

(δμa*,δμa*)=R̃(μa̅,μ̅s)[R(μa̅,μ̅s)R̃(μa̅,μ̅s)]1R(μa̅,μ̅s).

Here, R̃′(μ̅a,μ̅s) denotes the adjoint of R′(μ̅a, μ̅s). We prove in Appendix A that this minimum norm solution is a descent direction of the functional K(μa, μs).

Included in the minimum norm solution Eq. (9) is the operator [R′(μ̅a,μ̅s)R̃′(μ̅a,μ̅s)]-1. This operator is time consuming to calculate in practice because it requires solving as many direct and adjoint problems as there are data points. Instead, we approximate this operator by the scalar factor β multiplied to the identity operator, so that

(δμa*,δμa*)βR̃(μa̅,μ̅s)R(μa̅,μ̅s).

This approximation is substantially faster to compute. Moreover, we show in Appendix A that this approximate solution is a descent direction of the functional K(μa, μs) as well. To compute the update Eq. (10), we compute first r = R(μ̅a,μ̅s). This residual defined by Eq. (6) requires the solution of the direct problem with parameters μ̅a and μ̅s. Next, we must compute R̃′(μ̅a,μ̅s)r. This operation is given in terms of Ĩ, the solution of the adjoint problem [5, 6] given by

c1tĨ+ŝ·Ĩ+μaĨμsĨ=0inSn1×Ω×[0,T],
Ĩ=ronΓout={(ŝ,rb,t)Sn1×Ω×[0,T],ν(rb)·ŝ>0},
Ĩt=T=0onSn1×Ω.

This form for the adjoint problem relies on 𝓛 being self-adjoint. We consider this matter in detail in Appendix B.

Solving Eq. (11) is called backtransport. Hence, the transport-backtransport method is comprised of computing the direct problem followed by the adjoint problem. Upon solution of the direct and adjoint problems, we compute the components of Eq. (10) and thus obtain updated values given by

δμa*=β0TSn1Ĩ(ŝ,r,t)I(ŝ,r,t)dŝdt,
δμs*=β0TSn1Ĩ(ŝ,r,t)I(ŝ,r,t)dŝdt.

2.4. Procedure to solve the inverse problem

We now describe the procedure we take to solve the inverse problem. Suppose we have chosen an initial distribution for the absorption and scattering coefficients. Typically, we choose constant values for the initial distributions of μa and μs. Next, we choose an initial step size λ. Included within λ is the constant β given in Eq. (10). In the procedure below, this parameter λ changes to accelerate the convergence.

Once we have set these initial values, we proceed to solve the inverse problem in the following way.

  1. For each source in the set of sources, solve the direct problem comprised of Eq. (1) subject to Eq. (3) and Eq. (4) using the current distributions for μa and μs.
  2. Use the solution of the direct problem to compute M μa,μs (r d, t) defined by Eq. (5) for all of the detector locations r d ∈ ∂Ω and for t ∈ [0,T].
  3. Evaluate the cost functional K(μa,μs) given by Eq. (7).
    1. If this is the first iteration go to Step 4.
    2. If the difference of this value of K and the previous one is smaller than some user-defined tolerance, stop this procedure.
    3. If this value of K is smaller than the previous one, set λ ← 1.5λ and then go to Step 4.
    4. If this value of K is larger than the previous one, set λ ← 0.5λ, recalculate the parameter distributions through evaluation of Eq. (12) and return to Step 1.
  4. Using the residuals from each of the detectors, solve the adjoint problem Eq. (11) for Ĩ for each of the sources.
  5. For each source in the set of sources, calculate the updates for δμa and δμs using equation Eq. (12).
  6. Sum up all of these updates to obtain the final updates for δμa and δμs.
  7. Calculate the new parameter distributions by setting μaμa + λδμa and μsμs + λδμs.
  8. Return to Step 1.

It is possible that the line search method described above may affect the overall performance of a particular light scattering model in a way that is different for another light scattering model. Hence, using this algorithm in making absolute comparisons between the different light scattering models may be subject to objection. For that reason, we have taken care in analyzing the results that follow to avoid over-statement of the results.

3. Models for scattering operators

For the radiative transport equation Eq. (1), we must specify the scattering operator 𝓛. Tissues typically scatter light with a sharp forward-peak. For that case, the scattering phase function f given in Eq. (2) is is peaked about ŝ · ŝ ~ 1. Forward-peaked scattering makes the radiative transport equation difficult to solve. For that reason, we discuss several approximations that model forward-peaked scattering that we will consider in this study. These approximations are derived in detail in Ref. [4]. Here, we just present the approximations.

3.1. Henyey-Greenstein scattering

The Henyey-Greenstein scattering phase function,

fHG(ŝ·ŝ)=14π1g2(1+g22gŝ·ŝ)3/2=n=02n+14πgnPn(ŝ·ŝ),

is used typically for light propagation in tissues. Here, g denotes the asymmetry parameter or anisotropy factor and Pn denotes the nth order Legendre polynomial. For forward-peaked scattering, the anisotropy factor is g ~ 1. Using Eq. (13) in Eq. (2), we obtain 𝓛HG which we call the Henyey-Greenstein scattering operator.

3.2. The delta-Eddington approximation

The delta-Eddington approximation was introduced by Joseph et al. [16] as an extension to the two-term Eddington approximation to study forward-peaked and large-angle scattering. It approximates the scattering phase function as the sum of a delta function and a smooth background. When applying the delta-Eddington approximation to the Henyey-Greenstein scattering phase function, we obtain the following delta-Eddington scattering operator:

δEI=(1g2)I+1g24πs2[P0(ŝ·ŝ)+3g1+gP1(ŝ·ŝ)]I(ŝ,r)dŝ.

3.3. The Fokker-Planck approximation

The Fokker-Planck equation was first proposed by A. D. Fokker [17] and M. Planck [18] to describe the Brownian motion of a particle. In the context of light propagation in tissues, the Fokker-Planck approximation is given by [19]

FPI=12(1g)ΔŝI,

with ∆ŝ denoting the spherical laplacian.

3.4. The Fokker-Planck-Eddington approximation

Recently, the authors have combined the delta-Eddington approximation and the Fokker-Planck approximation to obtain the Fokker-Planck-Eddington approximation [4]:

FPEI=(1a0)I+a1ΔŝI+14πS2[b0P0(ŝ·ŝ)+3b1P1(ŝ·ŝ)]I(ŝr,t)dŝ.

For the Henyey-Greenstein scattering phase function defined in Eq. (13), the coefficients in Eq. (16) are given by

a0=g2(2g),
a1=g26(1g),
b0=12g2+g3,
b1=g3(35g+2g2).

3.5. The generalized Fokker-Planck-Eddington approximation

In addition to the Fokker-Planck-Eddington approximation, the authors have introduced the generalized Fokker-Planck-Eddington approximation [4]. This approximation uses a higher-order Fokker-Planck approximation along with the two-term Eddington approximation. This approximation leads to the generalized Fokker-Planck-Eddington scattering operator 𝓛gFPE defined as

gFPEI=(1a0)I+a1Δŝ(𝓘a2Δŝ)1I
+14πS2[b0P0(ŝ·ŝ)+3b1P1(ŝ·ŝ)]I(s'̂,r,t)dŝ',

with ℐ denoting the identity operator. For the Henyey-Greenstein scattering phase function defined in Eq. (13), the coefficients in Eq. (18) are given by

a0=g32g55g2
a1=7g36g2125g220g+4
a2=3g460g24
b0=2g45g3+5g25g2
b1=g8g327g2+27g827g8.

3.6. Two dimensional scattering

The approximations explained above are for the three dimensional scattering. For our computations, we consider only two spatial dimensions x and y and scattering is restricted to the xy-plane. For this case, the scattering phase function is slightly different and so are the corresponding approximations.

According to J. Heino et al [20], the Henyey-Greenstein phase function in two dimensions is:

fHG=fHG(θθ)=12π1g21+g22gcos(θθ)=12π+1πn=1gncos[n(θθ)].

The corresponding Henyey-Greenstein scattering operator is given by

HGI=I+02πfHG(θθ)I(θ',x,y,t)dθ.

Consequently, the approximations for the scattering operator are given by the following.

  • delta-Eddington scattering operator

    δEI=(1g2)I+1g22π02π[1+2g1+gcos(θθ)]I(θ',x,y,t)dθ.

  • Fokker-Planck scattering operator

    FPI=(1g)θ2I.

  • Fokker-Planck-Eddington scattering operator

    FPEI=(1a0)I+a1θ2I+12π02π[b0+b1cos(θθ)]I(θ',x,y,t)dθ,

    with

    a0=g25(94g),

    a1=g25(1g),

    b0=15(59g2+4g3),

    b1=g5(58g+3g2).

  • generalized Fokker-Planck-Eddington scattering operator

    gFPEI=(1a0)I+a1θ2(Ia2θ2)1I
    +12π02π[b0+b1cos(θθ)]I(θ,x,y,t)dθ,

    with

    a0=g3207g20g7,

    a1=105g31g21600g21120g+196,

    a2=75g80g28,

    b0=7g420g3+20g720g7,

    b1=g98g4389g3+501g2−259g+49340g2259g+49.

4. Numerical solution of the direct scattering problem

We describe the numerical solution of

c1tI+cosθxI+sinθyI+μaIμsI=Q,

in [0,2π]×[0,Lx]×[0,Ly]×[0,T]. Here, we describe the numerical solution for a general scattering operator 𝓛. Later, we give numerical approximations that we use for computing numerical approximations of 𝓛HG defined by Eq. (21), 𝓛δE defined by Eq. (22), 𝓛FP defined by Eq. (23), 𝓛FPE defined by Eq. (24) and 𝓛gFPE defined by Eq. (26).

4.1. Source iteration method

To solve Eq. (28), we use the source iteration method [21]. For the source iteration method, we solve numerically

c1tIl+cosθxIl+sinθyIl+μaIl=Jl1

for l = 1,2,⋯ until we reach convergence:

IlIl1<ε,

with ε denoting a user-defined termination parameter. Here, the solution we seek I is given by Il as l → ∞. The source Jl is defined as

J0=Ql=0,
Jl=μsIl+Q,l=1,2,.

For each iteration, the numerical method is comprised of two parts: computing the source update Eq. (31) and solving the partial differential equation Eq. (29). In what follows, we describe each of these two parts in detail.

4.2. Computing the source update

Depending on the scattering operator we are considering, the source update involves either an integral operator, a differential operator or a combination of both. For integral operators, we use a quadrature rule and for the differential operators, we use a finite difference scheme.

For both types of operators, we consider the angular grid of N points defined as θn = (n - 1)∆θ with ∆θ = 2π/N for n = 1, ⋯,N. We now introduce the approximant In(x,y,t) ≈ I(θn,x,y,t) that we will use in the numerical approximations below.

For integral operations, we use a product quadrature rule [22] of the form:

02πf(θnθ)I(θ,x,y,t)dθk=1NfnkIk(x,y,t),n=1,,N.

Here, the weights fnk are defined as

fnk=θkΔθ/2θk+Δθ/2f(θnθ)dθ.

For the derivative operations, we use a centered finite difference scheme of the form

θ2InIn12In+In+1(Δθ)2,n=1,,N,

with periodic numerical boundary conditions: I -1 = IN and I N+1 = I 1.

By using either Eq. (31), Eq. (33) or a combination of the two, we form the discrete scattering operator 𝓛N. Using 𝓛N, we compute the source update at θn through evaluation of

Jn=μs𝓛NIn(x,y,t)+Q(θn,x,y,t).

4.3. Diamond differencing scheme

Using the angle grid defined above, we seek to solve numerically the partial differential equation

c1tIn+cosθnxIn+sinθnyIn+μaIn=Jn,n=1,,N,

with Jn known. To solve this problem, we use the diamond differencing scheme [21]. An equispaced spatial grid defines spatial cells each of size ∆x by ∆y. We use integer indices i and j to denote each cell. Furthermore, we assume that μa and μs are constant within each cell. Then, we seek the cell-centered average of In(x,y,t) over each cell defined as

Inij(t)=1ΔxΔyxi1/2xi+1/2yi1/2yi+1/2In(x,y,t)dydx.

Here, we have denoted the edges of spatial cells by the fractional indices i ± 1/2 and j ± 1/2.

To derive an equation for Inij(t), we integrate Eq. (36) over [x i-1/2,x i+1/2] × [y j-1/2,y j+1/2], divide by ∆xy and obtain the semi-discrete system

c1tInij+1Δxcosθn(In,i+1/2,jIn,i1/2,j)+1Δysinθn(In,i,j+1/2In,i,j1/2)
+μa,ijInij=Jnij.

Here, we have introduced the cell-edge quantities:

In,i±1/2,j=1Δyyj1/2yj+1/2In(xi±1/2,y,t)dy,
In,i,j±1/2=1Δxxi1/2xi+1/2In(x,yj±1/2,t)dx.

Next, we seek a method to evolve the solution from time tk to t k+1 = tk+∆t. Integrating Eq. (38) with respect to t over [t k-1/2,t k+1/2] and dividing by ∆t, we obtain

1cΔt(Inijk+1/2Inijk1/2)+1Δxcosθn(In,i+1/2,jkIn,i1/2,jk)
+1Δysinθn(In,i,j+1/2kIn,i,j1/2k)+μa,ijInijk=Jnijk.

Here, we have introduced the notation Inij k+1/2Inij(t k+1/2), Inij k-1/2Inij(t k-1/2) and

Inijk=1Δttk1/2tk+1/2Inij(t)dt.

The diamond differencing scheme uses the following approximations to go back-and-forth from cell-averaged quantities (fractional indices) to cell-edge quantities (integer indices):

Inijk=12(In,i+1/2,jk+In,i1/2,jk),
Inijk=12(In,i,j+1/2k+In,i,j1/2k),
Inijk=12(Inijk+1/2+Inijk1/2).

By combining Eq. (40) and Eq. (42), we obtain a discrete system of equations whose solution yields the numerical solution we seek. For the diamond differencing scheme in time, we substitute 2Inijk - Inij k+1/2 for Inij k-1/2 in Eq. (40) and obtain an update formula for Inij k+1/2. Once we determine that update, we compute the full update through evaluation of Inij k+1 = 2Inij k+1/2 -Inij k. The actual form of the update formula for Inij k+1/2 depends on the signs of cos θn and sin θn. In particular, we use the following.

  • For cos θn ≥ 0 and sin θn ≥ 0, we make the substitutions

    I n,i+1/2,j → 2Inij - I n,i-1/2,j and I n,i,j+1/2 → 2Inij - I n,i,j-1/2.

  • For cos θn ≥ 0 and sin θn ≤ 0, we make the substitutions

    I n,i+1/2,j ← 2Inij - I n,i-1/2,j and I n,i,j-1/2 ← 2Inij - I n,i,j+1/2.

  • For cos θn ≤ 0 and sin θn ≥ 0, we make the substitutions

    I n,i-1/2,j ← 2Inij - I n,i+1/2,j and I n,i,j+1/2 ← 2Inij - I n,i,j-1/2.

  • For cos θn ≤ 0 and sin θn ≤ 0, we make the substitutions

    I n,i-1/2,j ← 2Inij - I n,i+1/2,j and I n,i,j+1/2 ← 2Inij - I n,i,j+1/2.

So for example, we solve the following system of equations to compute the Inij k+1/2 update for cos θn ≥ 0 and sin θn ≥ 0:

Inijk+1/2=InijkcΔtΔxcosθn(InijkIn,i1/2,jk)cΔtΔysinθn(InijkIn,i,j1/2k)
cΔt2μa,ijInijk+cΔt2Jnijk.

To ensure stable time-stepping, we must set cΔt12min{x,y}.. In addition, we must set ctmaxij{μa,ij} < 1 [5, 21]. To avoid this second restriction, we have implemented a semi-implicit time-step. Thus, instead of solving Eq. (43), we solve

(1+cΔt2μa,ij)Inijk+1/2=InijkcΔtΔxcosθn(InijkIn,i1/2,jk)
cΔtΔysinθn(InijkIn,i,j1/2k)+cΔt2Jnijk.

5. Numerical results

We investigate the quality of image reconstructions produced by the five light scattering models: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation. For all numerical experiments, we obtain simulated boundary measurements by solving the radiative transport equation with the Henyey-Greenstein scattering phase function with g = 0.9. Using those simulated data, we use the transport-backtransport method to reconstruct images of the absorption and scattering coefficients. We perform the transport-backtransport method using each of the five light scattering models to obtain five sets of image reconstructions.

For our numerical calculations, the spatial domain is a 2 cm × 2 cm square with constant background optical properties: μa = 0.1mm-1 and μs = 1.0 mm-1. Within this medium, there is an absorbing obstacle whose absorption coefficient is 1.0 mm-1 and a scattering obstacle whose scattering coefficient is 5.0mm-1. These parameter values are reasonable for studies of light propagation in tissues. In particular, we have chosen these parameters to focus on the regime in which the scattering mean free path ls = μs -1 is significantly smaller than the transport mean free path I * = [μs(1 - g)]-1. For that case, the medium scatters light strongly, but with a strong forward peak. This regime is where the approximations are designed to work well. Figures 1(a) and 2(a) show the placement and size of the absorbing and scattering obstacles, respectively.

For all of our computations of the direct and adjoint problems, we have set the number of discrete angles to be N = 64 so that ∆θ = π/32. The spatial grid has 51 points in x and y so that ∆x = ∆y = 0.04cm. We have taken a time step of ∆t = 0.1 sampling an 8 second interval of time (we have set c = 1 here for convenience).

Along the boundary of the spatial domain, we have placed forty, equi-distant sources and detectors surrounding the entire spatial domain. In fact, these sources and detectors are collocated spatially. Thus, for a particular source, the set of detectors correspond to all source/detector locations except the ones located in the same side of the domain.

The results of these numerical simulations are shown in Figs. 1 and 2. In all of the reconstructed images, we observe the presence of “cross-talk” in which the scattering obstacle appears as an artifact in the absorption coefficient image and the absorbing obstacle appears as an artifact in the scattering coefficient image. Even the so-called inverse crime solution in which the radiative transport equation is used as the light scattering model exhibits cross-talk. This cross-talk is due to the ill-posedness of the inverse problem. In particular, the data we have are limited to just a few spatial locations (albeit surrounding the entire domain) of angle-averaged quantities. Using a priori information or incorporating more data into the inverse problem may alleviate this problem.

Even though all of the light scattering models produce images with substantial cross-talk, all of the reconstructed images are reasonably good. For example, the absorbing and scattering obstacles are localized correctly. We observe that the reconstructions of the scattering coefficient are better overall than the respective images for the absorption coefficient. Differences occurring between the images reconstructed using the different light scattering models include better localization and less image artifacts. For example, the reconstructions for delta-Eddington, Fokker-Planck and Fokker-Planck-Eddington models show image artifacts near the source positions. Additional regularization near the sources may remove these artifacts provided we the reasonable assumption that the obstacles to be reconstructed are not close to the sources.

For the absorption coefficient image reconstructions shown in Fig. 1, the cross-talk is more severe than for the scattering coefficient. Not surprisingly, the inverse crime solution in which the image is produced by the radiative transport equation is the best one. The worst qualitative and quantitative image is produced by the Fokker-Planck approximation. The delta-Eddington approximation produces a slightly better image and the Fokker-Planck-Eddington approximation produces an even better image. Besides that inverse crime solution, the generalized Fokker-Planck-Eddington scattering model produced the best image.

For the scattering coefficient, we find again that the inverse crime solution produces the best image. Besides that inverse crime solution, we find that the generalized-Fokker-Planck-Eddington approximation produces the best image. The Fokker-Planck approximation produces a slightly better qualitative result than the delta-Eddington approximation since there are less image artifacts. However, the delta-Eddington approximation yields a better quantitative result. The Fokker-Planck-Eddington approximation is better than delta-Eddington and Fokker-Planck, but worse than the generalized Fokker-Planck-Eddington approximation.

In Fig. 5 we show the evolution of the cost functional Eq. (7) as a function of the successful iterations using the different light scattering models. This plot gives another perspective on the performance of the different light scattering models. For example, we see clearly that the radiative transport equation performs best by a substantial amount. The wide variety of cost functional values are due largely to modeling error. When the light scattering model agrees better with the radiative transport equation, we find that the cost values are smaller than those from light scattering models that do not agree as well. That is the reason why we find that the generalized Fokker-Planck-Eddington cost values are smaller than the delta-Eddington or Fokker-Planck ones. In fact, we find that the generalized Fokker-Planck-Eddington approximation performs best among the approximations. In these computations, we have set a limit on the number of consecutive iterations in which the residual does not decrease rather than setting a tolerance value (Step 3 (b) in the procedure given in Section 2.4). It is for that reason that we see that the delta-Eddington approximation gives the poorest performance in this context – it reaches the maximum number of rejected iterations faster than the others without reducing the cost functional substantially.

 

Fig. 1. The actual absorption coefficient (a) and reconstructions of the absorption coefficient produced by the five light scattering models (b)–(f).

Download Full Size | PPT Slide | PDF

 

Fig. 2. The actual scattering coefficient (a) and reconstructions of the scattering coefficient produced by the five light scattering models (b)–(f).

Download Full Size | PPT Slide | PDF

 

Fig. 3. Comparison of the values of the cost functional as a function of successful iterations for the radiative transport equation (RTE), the generalized Fokker-Planck-Eddington approximation (gFPE), the Fokker-Planck-Eddington approximation (FPE), the Fokker-Planck approximation (FP) and delta-Eddington approximation (DE).

Download Full Size | PPT Slide | PDF

Tables Icon

Table 1. Time required to compute the numerical solution of a single direct problem in seconds.

An important consideration in using one of the various approximations for image reconstruction is the computational savings they afford over using the radiative transport equation. For the computational results shown here, we have used a workstation with two dual core Intel Xeon (3.0 GHz) processors. In Table 1 we give the time required to solve a single direct problem for each of the light scattering models. Using any of the approximations instead of the radiative transport equation is at least 25% faster. This computational savings is compounded even further since the transport-backtransport method requires the computation of several direct and adjoint problems. The implementation of the numerical solution may be implemented on parallel cluster thereby making these calculations faster overall. For example, one may run each solution corresponding to an individual source on different processors. Doing so would nearly linearly with the number of processors. However, the relative speeds among the different light scattering models are still an important factor in the overall computational time.

6. Conclusions

We have investigated the use of approximations of the radiative transport equation with forward-peaked scattering for producing image reconstructions. Using the transport-backtransport method for solving the inverse problem, we have used five different light scattering models: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation to produce five sets of image reconstructions. By comparing those results, we have evaluated the potential use of approximations for reconstructing images in diffuse optical tomography. These approximations provide substantial computational savings compared with the computational resources required by the radiative transport equation.

All of our results exhibit cross-talk in which the absorbing obstacle appears in the scattering coefficient image and the scattering obstacle appears in the absorption coefficient image. Crosstalk occurs even in the inverse crime scenario in which the radiative transport equation computes the data and is the light scattering model for the inverse problem solution. It is an artifact of the ill-posedness inherent in the inverse problem. Hence, it has become clear through this investigation that eliminating or reducing cross-talk, regardless of the light scattering model used in the inverse solution, is an important topic. We will address cross-talk in our future work.

Overall, we have found that all of the light scattering models produce reasonable results. Besides the inverse crime solution that uses the radiative transport equation, we have found that the generalized Fokker-Planck-Eddington approximation produces the best qualitative and quantitative images. As expected, the quality of image reconstructions is directly proportional to the computational time required to solve a single direct problem. Nonetheless, since the approximations are at least 25% faster to solve than the radiative transport equation, they are certainly worthwhile to consider for practical applications.

Appendix A

Although the proof that Eq. (9) and Eq. (10) are in the descent direction of K given by Eq. (7) can be found in Refs. [13] and [14], we give the proof here for reference. Using the definition of K we have that

K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R((μ̅a,μ̅s)+ω(δμa*,δμs*)L22.

Using Eq. (8) we rewrite Eq. (A1) as

K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R(μ̅a,μ̅s)+R(μ̅a,μ̅s)[ω(δμa*,δμs*)]
+O((δμa*,δμs*)L22)L22.

By using the definition of the L 2 norm, we find that

K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R(μ̅a,μ̅s),R(μ̅a,μ̅s)
+R(μ̅a,μ̅s),R(μ̅a,μ̅s)[ω(δμa*,δμs*)]+O((δμa*,δμs*)L22)

Now, using the definition of (δμa *,δμs *) given by Eq. (9), we obtain

K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R(μ̅a,μ̅s),R(μ̅a,μ̅s)
+R(μ̅a,μ̅s),ωR(μ̅a,μ̅s)R̃(μ̅a,μ̅s)[R(μ̅a,μ̅s)R̃(μ̅a,μ̅s)]1R(μ̅a,μ̅s)
+O((δμa*,δμs*)L22).

Simplifying the second term on the right hand side of Eq. (A4) and applying the definition of the L 2 norm again, we finally arrive at

K((μ̅a,μ̅s)+ω(δμa*,δμs*))=K(μ̅a,μ̅s)ωR(μ̅a,μ̅s)L22+O(δμa*,δμs*L22).

Since ω is a positive constant, we have thus proved that (δμa *, δμs *) is a descent direction of K(μa,μs).

We can also prove that the approximation given by Eq. (10) is also a descent direction of K(μa,μs). By substituting -β R̃′(μ̅a, μ̅s)R(μ̅a, μ̅s) in for (δμa *, δμs *) in Eq. (A3), we obtain

K((μ̅a,μ̅s)+ω(δμa*,δμs*))=K(μ̅a,μ̅s)ωβR'̃(μ̅a,μ̅s)R(μ̅a,μ̅s)L22
+O(δμa*,δμs*L22),

with ωβ a positive constant. Thus, we have proved also that Eq. (10) is a descent direction of K(μa,μs).

Appendix B

The adjoint problem given in Eq. (11) relies on the scattering operator 𝓛 being self-adjoint. By self-adjoint, we mean that 𝓛 satisfies

Sn1Ĩ(ŝ)[I(ŝ)]dŝ=Sn1[Ĩ(ŝ)]I(ŝ)dŝ,

for any two suitable functions I and Ĩ defined on S n-1. For example, the identity operator ℐ is self-adjoint.

For any scattering phase function that depends only on ŝ · ŝ′, the corresponding scattering operator is self-adjoint. Thus, the Henyey-Greenstein (Section 3.1) and delta-Eddington (Section 3.2) scattering operators are self-adjoint. Moreover, the integral terms in the Fokker-Planck-Eddington (Section 3.4) and the generalized Fokker-Planck-Eddington (Section 3.5) scattering operators are self-adjoint. By a straight-forward application of partial integration, one can show readily that the spherical laplacian ∆ŝ is self-adjoint.

It remains to show that 𝓛gFP I = a 1ŝ(ℐ - a 2ŝ)-1 I is self-adjoint. To do so, consider the inner product

Ĩ,gFPI=Sn1Ĩ(ŝ)[a1Δŝ(a2Δŝ)1I(ŝ)]dŝ.

Let I(ŝ)=(ℐ-a 2ŝ)J(ŝ) and Ĩ(ŝ)=(ℐ-a 2s ̂)J̃(ŝ). Then,

Ĩ,gFPI=Sn1[(a2Δŝ)J̃(ŝ)][a1ΔŝJ(ŝ)]dŝ
=Sn1[a1ΔŝJ̃(ŝ)][(a2Δŝ)J(ŝ)]dŝ
=Sn1[a1Δŝ(a2Δŝ)1Ĩ(ŝ)]I(ŝ)dŝ
=gFPI,Ĩ.

Thus, 𝓛gFP is self-adjoint. Since all of the components of the light scattering models are self-adjoint, all of the approximations are self-adjoint. Therefore, the adjoint problem has the form given in Eq. (11).

Acknowledgments

P. González-Rodríguez and A. D. Kim acknowledge support from the National Science Foundation (grants DMS-0504858 and EEC-0616228) and the US Department of Energy through the University of California, Merced Center for Computational Biology. P. González-Rodríguez also acknowledges support from the Spanish Ministry of Education and Science (Project FIS2007-62673).

References and links

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

2. A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005). [CrossRef]  

3. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, 1997).

4. P. Gonzàlez-Rodríguez and A. D. Kim,, “Light propagation in tissues with forward-peaked and large-angle scattering,” Appl. Opt. 47, 2599–2609 (2008). [CrossRef]   [PubMed]  

5. O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998). [CrossRef]  

6. O. Dorn, “Scattering and absorption transport sensitivity functions for optical tomography,” Opt. Express 13, 492–506 (2000). [CrossRef]  

7. F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Probl. 11, 1225–1232 (1995). [CrossRef]  

8. M. Vogeler, “Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation-backpropagation method,” Inverse Probl. 19, 739–753 (2003). [CrossRef]  

9. O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999). [CrossRef]  

10. O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000). [CrossRef]  

11. P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005). [CrossRef]  

12. E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000). [CrossRef]  

13. F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, PA, 2001). [CrossRef]  

14. C. R. Vogel, Computational Methods for Inverse Problems (SIAM, Philadelphia, PA, 2002). [CrossRef]  

15. F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002). [CrossRef]  

16. J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976). [CrossRef]  

17. A. D. Fokker, “Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld,” Ann. Phys. 43, 810–820 (1914). [CrossRef]  

18. M. Planck, “Ueber einen Satz der statistichen Dynamik und eine Erweiterung in der Quantumtheorie,” Sitzung-ber. Preuss. Akad. Wiss. 5, 324–341 (1917).

19. A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20, 92–98 (2003). [CrossRef]  

20. J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003). [CrossRef]  

21. E. E. Lewis and W. F. Miller Jr., Computational Methods of Neutron Transport (American Nuclear Society, Inc., La Grange Park, IL, 1993).

22. L. M. Delves and J. L. Mohamed, Computational Methods for Integral Equations (Cambridge University Press, Cambridge, 1985). [CrossRef]  

References

  • View by:
  • |
  • |
  • |

  1. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).
    [Crossref]
  2. A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005).
    [Crossref]
  3. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, 1997).
  4. P. Gonzàlez-Rodríguez and A. D. Kim,, “Light propagation in tissues with forward-peaked and large-angle scattering,” Appl. Opt. 47, 2599–2609 (2008).
    [Crossref] [PubMed]
  5. O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998).
    [Crossref]
  6. O. Dorn, “Scattering and absorption transport sensitivity functions for optical tomography,” Opt. Express 13, 492–506 (2000).
    [Crossref]
  7. F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Probl. 11, 1225–1232 (1995).
    [Crossref]
  8. M. Vogeler, “Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation-backpropagation method,” Inverse Probl. 19, 739–753 (2003).
    [Crossref]
  9. O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
    [Crossref]
  10. O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000).
    [Crossref]
  11. P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
    [Crossref]
  12. E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000).
    [Crossref]
  13. F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, PA, 2001).
    [Crossref]
  14. C. R. Vogel, Computational Methods for Inverse Problems (SIAM, Philadelphia, PA, 2002).
    [Crossref]
  15. F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
    [Crossref]
  16. J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976).
    [Crossref]
  17. A. D. Fokker, “Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld,” Ann. Phys. 43, 810–820 (1914).
    [Crossref]
  18. M. Planck, “Ueber einen Satz der statistichen Dynamik und eine Erweiterung in der Quantumtheorie,” Sitzung-ber. Preuss. Akad. Wiss. 5, 324–341 (1917).
  19. A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20, 92–98 (2003).
    [Crossref]
  20. J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
    [Crossref]
  21. E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport (American Nuclear Society, Inc., La Grange Park, IL, 1993).
  22. L. M. Delves and J. L. Mohamed, Computational Methods for Integral Equations (Cambridge University Press, Cambridge, 1985).
    [Crossref]

2008 (1)

2005 (2)

P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
[Crossref]

A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005).
[Crossref]

2003 (3)

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
[Crossref]

A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20, 92–98 (2003).
[Crossref]

M. Vogeler, “Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation-backpropagation method,” Inverse Probl. 19, 739–753 (2003).
[Crossref]

2002 (1)

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

2000 (3)

O. Dorn, “Scattering and absorption transport sensitivity functions for optical tomography,” Opt. Express 13, 492–506 (2000).
[Crossref]

E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000).
[Crossref]

O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000).
[Crossref]

1999 (2)

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

O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
[Crossref]

1998 (1)

O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998).
[Crossref]

1995 (1)

F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Probl. 11, 1225–1232 (1995).
[Crossref]

1976 (1)

J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976).
[Crossref]

1917 (1)

M. Planck, “Ueber einen Satz der statistichen Dynamik und eine Erweiterung in der Quantumtheorie,” Sitzung-ber. Preuss. Akad. Wiss. 5, 324–341 (1917).

1914 (1)

A. D. Fokker, “Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld,” Ann. Phys. 43, 810–820 (1914).
[Crossref]

Arridge, S.

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
[Crossref]

Arridge, S. R.

A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005).
[Crossref]

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

Ascher, U.

E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000).
[Crossref]

Berryman, J. G.

O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
[Crossref]

Bertete-Aguirre, H.

O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
[Crossref]

Delves, L. M.

L. M. Delves and J. L. Mohamed, Computational Methods for Integral Equations (Cambridge University Press, Cambridge, 1985).
[Crossref]

Dierkes, T.

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

Dorn, O.

P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
[Crossref]

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000).
[Crossref]

O. Dorn, “Scattering and absorption transport sensitivity functions for optical tomography,” Opt. Express 13, 492–506 (2000).
[Crossref]

O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
[Crossref]

O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998).
[Crossref]

Fokker, A. D.

A. D. Fokker, “Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld,” Ann. Phys. 43, 810–820 (1914).
[Crossref]

Gibson, A. P.

A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005).
[Crossref]

Gonzàlez-Rodríguez, P.

P. Gonzàlez-Rodríguez and A. D. Kim,, “Light propagation in tissues with forward-peaked and large-angle scattering,” Appl. Opt. 47, 2599–2609 (2008).
[Crossref] [PubMed]

P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
[Crossref]

Haber, E.

E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000).
[Crossref]

Hebden, J. C.

A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005).
[Crossref]

Heino, J.

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
[Crossref]

Ishimaru, A.

A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, 1997).

Joseph, J. H.

J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976).
[Crossref]

Keller, J. B.

Kim, A. D.

Kim,, A. D.

Kindelan, M.

P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
[Crossref]

Lewis, E. E.

E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport (American Nuclear Society, Inc., La Grange Park, IL, 1993).

Miller, E. L.

O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000).
[Crossref]

Miller, W. F.

E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport (American Nuclear Society, Inc., La Grange Park, IL, 1993).

Mohamed, J. L.

L. M. Delves and J. L. Mohamed, Computational Methods for Integral Equations (Cambridge University Press, Cambridge, 1985).
[Crossref]

Moscoso, M.

P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
[Crossref]

Natterer, F.

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Probl. 11, 1225–1232 (1995).
[Crossref]

F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, PA, 2001).
[Crossref]

Oldenburg, D.

E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000).
[Crossref]

Palamodov, V.

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

Papanicolaou, G. C.

O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
[Crossref]

Planck, M.

M. Planck, “Ueber einen Satz der statistichen Dynamik und eine Erweiterung in der Quantumtheorie,” Sitzung-ber. Preuss. Akad. Wiss. 5, 324–341 (1917).

Rappaport, C.

O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000).
[Crossref]

Sielschott, H.

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

Sikora, J.

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
[Crossref]

Somersalo, E.

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
[Crossref]

Vogel, C. R.

C. R. Vogel, Computational Methods for Inverse Problems (SIAM, Philadelphia, PA, 2002).
[Crossref]

Vogeler, M.

M. Vogeler, “Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation-backpropagation method,” Inverse Probl. 19, 739–753 (2003).
[Crossref]

Weinman, J. A.

J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976).
[Crossref]

Wiscombe, W. J.

J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976).
[Crossref]

Wübbeling, F.

F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Probl. 11, 1225–1232 (1995).
[Crossref]

F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, PA, 2001).
[Crossref]

Ann. Phys. (1)

A. D. Fokker, “Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld,” Ann. Phys. 43, 810–820 (1914).
[Crossref]

Appl. Opt. (1)

Inverse Probl. (8)

O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998).
[Crossref]

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

F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Probl. 11, 1225–1232 (1995).
[Crossref]

M. Vogeler, “Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation-backpropagation method,” Inverse Probl. 19, 739–753 (2003).
[Crossref]

O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl. 15, 1523–1558 (1999).
[Crossref]

O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adoint fields and level sets,” Inverse Probl. 16, 1119–1156 (2000).
[Crossref]

P. Gonzàlez-Rodríguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. 21, 565–590 (2005).
[Crossref]

E. Haber, U. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” Inverse Probl. 16, 1263–1280 (2000).
[Crossref]

J. Atmos. Sci. (1)

J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. 33, 2452–2459 (1976).
[Crossref]

J. Opt. Soc. Am. A (1)

Opt. Express (1)

O. Dorn, “Scattering and absorption transport sensitivity functions for optical tomography,” Opt. Express 13, 492–506 (2000).
[Crossref]

Phys. Med. Bio. (1)

A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Bio. 50, R1–R43 (2005).
[Crossref]

Phys. Rev. E (1)

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003).
[Crossref]

SIAM J. Appl. Math. (1)

F. Natterer, H. Sielschott, O. Dorn, T. Dierkes, and V. Palamodov, “Fréchet derivatives for some bilinear inverse problems,” SIAM J. Appl. Math. 62, 2092–2113 (2002).
[Crossref]

Sitzung-ber. Preuss. Akad. Wiss. (1)

M. Planck, “Ueber einen Satz der statistichen Dynamik und eine Erweiterung in der Quantumtheorie,” Sitzung-ber. Preuss. Akad. Wiss. 5, 324–341 (1917).

Other (5)

A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, 1997).

F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, PA, 2001).
[Crossref]

C. R. Vogel, Computational Methods for Inverse Problems (SIAM, Philadelphia, PA, 2002).
[Crossref]

E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport (American Nuclear Society, Inc., La Grange Park, IL, 1993).

L. M. Delves and J. L. Mohamed, Computational Methods for Integral Equations (Cambridge University Press, Cambridge, 1985).
[Crossref]

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (3)

Fig. 1.
Fig. 1. The actual absorption coefficient (a) and reconstructions of the absorption coefficient produced by the five light scattering models (b)–(f).
Fig. 2.
Fig. 2. The actual scattering coefficient (a) and reconstructions of the scattering coefficient produced by the five light scattering models (b)–(f).
Fig. 3.
Fig. 3. Comparison of the values of the cost functional as a function of successful iterations for the radiative transport equation (RTE), the generalized Fokker-Planck-Eddington approximation (gFPE), the Fokker-Planck-Eddington approximation (FPE), the Fokker-Planck approximation (FP) and delta-Eddington approximation (DE).

Tables (1)

Tables Icon

Table 1. Time required to compute the numerical solution of a single direct problem in seconds.

Equations (88)

Equations on this page are rendered with MathJax. Learn more.

c1tI+ŝ·I+μaIμsI=QonSn1×Ω×[0,T].
I=I+Sn1f(ŝ·ŝ)I(ŝ,r,t)dŝ.
I=0onΓin={(ŝ,rb,t)Sn1×Ω×[0,T],ν(rb)·ŝ<0},
It=0 =0onSn1×Ω.
Mμa,μs (rb,t) = Soutn1ŝ.ν(rd)I(ŝ,rd,t)dŝonΩ×[0,T].
R:P×PD,R(μa,μs)=Mμa,μsF.
K(μa,μs)=12R(μa,μs)L22.
R (μa,μs)R(μa̅,μ̅s)+R'(μa̅,μ̅s)(δμa,δμs),
(δμa*,δμa*)=R̃(μa̅,μ̅s)[R(μa̅,μ̅s)R̃(μa̅,μ̅s)]1R(μa̅,μ̅s).
(δμa*,δμa*)βR̃(μa̅,μ̅s)R(μa̅,μ̅s).
c1tĨ+ŝ·Ĩ+μaĨμsĨ=0inSn1×Ω×[0,T],
Ĩ=ronΓout={(ŝ,rb,t)Sn1×Ω×[0,T],ν(rb)·ŝ>0},
Ĩt=T =0onSn1×Ω.
δμa*=β0TSn1Ĩ(ŝ,r,t)I(ŝ,r,t)dŝdt,
δμs*=β0TSn1Ĩ(ŝ,r,t)I(ŝ,r,t)dŝdt.
fHG(ŝ·ŝ)=14π1g2(1+g22gŝ·ŝ)3/2=n=02n+14πgnPn(ŝ·ŝ),
δEI=(1g2)I+1g24πs2[P0(ŝ·ŝ)+3g1+gP1(ŝ·ŝ)]I(ŝ,r)dŝ.
FPI=12(1g)ΔŝI,
FPE I =(1a0)I+a1ΔŝI+14πS2[b0P0(ŝ·ŝ)+3b1P1(ŝ·ŝ)]I(ŝr,t)dŝ.
a0=g2(2g),
a1=g26(1g) ,
b0=12g2+g3,
b1=g3(35g+2g2).
gFPE I=(1a0)I+a1Δŝ(𝓘a2Δŝ)1I
+14π S2[b0P0(ŝ·ŝ)+3b1P1(ŝ·ŝ)]I(s'̂,r,t)dŝ',
a0=g3 2g55g2
a1=7g36 g2125g220g+4
a2=3g460g24
b0=2g45g3+5g25g2
b1=g8g327g2+27g827g8 .
fHG=fHG(θθ)=12π1g21+g22gcos(θθ)=12π+1πn=1gncos[n(θθ)].
HGI=I+02πfHG(θθ)I(θ',x,y,t)dθ.
δEI=(1g2)I+1g22π02π[1+2g1+gcos(θθ)]I(θ',x,y,t)dθ.
FPI=(1g)θ2I.
FPEI=(1a0)I+a1θ2I+12π02π[b0+b1cos(θθ)]I(θ',x,y,t)dθ,
a0=g25 (94g) ,
a1=g25(1g),
b0=15 (59g2+4g3) ,
b1=g5(58g+3g2).
gFPEI=(1a0)I+a1θ2(Ia2θ2)1I
+12π 02π[b0+b1cos(θθ)] I (θ,x,y,t) d θ ,
a0=g3207g20g7,
a1=105g31g21600g21120g+196,
a2=75g80g28,
b0=7g420g3+20g720g7,
b1=g98g4389g3+501g2−259g+49340g2259g+49.
c1tI+cosθxI+sinθyI+μaIμsI=Q,
c1tIl+cosθxIl+sinθyIl+μaIl=Jl1
IlIl1<ε ,
J0=Ql=0,
Jl=μsIl+Q,l=1,2,.
02πf(θnθ)I(θ,x,y,t)dθk=1NfnkIk(x,y,t),n=1,,N.
fnk=θkΔθ/2θk+Δθ/2f(θnθ)dθ.
θ2 InIn12In+In+1(Δθ)2,n=1,,N,
Jn=μs𝓛NIn(x,y,t)+Q(θn,x,y,t).
c1tIn+cosθnxIn+sinθnyIn+μaIn=Jn,n=1,,N,
Inij(t)=1ΔxΔyxi1/2xi+1/2yi1/2yi+1/2In(x,y,t)dydx.
c1tInij+1Δx cos θn (In,i+1/2,jIn,i1/2,j) +1Δysinθn(In,i,j+1/2In,i,j1/2)
+μa,ijInij=Jnij.
In,i±1/2,j=1Δyyj1/2yj+1/2In(xi±1/2,y,t)dy,
In,i,j±1/2=1Δxxi1/2xi+1/2In(x,yj±1/2,t)dx.
1cΔt(Inijk+1/2Inijk1/2)+1Δxcosθn(In,i+1/2,jkIn,i1/2,jk)
+1Δysinθn(In,i,j+1/2kIn,i,j1/2k)+μa,ijInijk=Jnijk.
Inijk=1Δttk1/2tk+1/2Inij(t)dt.
Inijk=12(In,i+1/2,jk+In,i1/2,jk),
Inijk=12(In,i,j+1/2k+In,i,j1/2k),
Inijk=12(Inijk+1/2+Inijk1/2).
Inijk+1/2=InijkcΔtΔxcosθn(InijkIn,i1/2,jk)cΔtΔysinθn(InijkIn,i,j1/2k)
cΔt2μa,ijInijk+cΔt2Jnijk.
(1+cΔt2μa,ij)Inijk+1/2=InijkcΔtΔxcosθn(InijkIn,i1/2,jk)
cΔtΔysinθn(InijkIn,i,j1/2k)+cΔt2Jnijk.
K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R((μ̅a,μ̅s)+ω(δμa*,δμs*)L22.
K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R(μ̅a,μ̅s)+R(μ̅a,μ̅s)[ω(δμa*,δμs*)]
+O((δμa*,δμs*)L22)L22.
K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R(μ̅a,μ̅s),R(μ̅a,μ̅s)
+R(μ̅a,μ̅s),R(μ̅a,μ̅s)[ω(δμa*,δμs*)]+O((δμa*,δμs*)L22)
K((μ̅a,μ̅s)+ω(δμa*,δμs*))=12R(μ̅a,μ̅s),R(μ̅a,μ̅s)
+R(μ̅a,μ̅s),ωR(μ̅a,μ̅s)R̃(μ̅a,μ̅s)[R(μ̅a,μ̅s)R̃(μ̅a,μ̅s)]1R(μ̅a,μ̅s)
+O((δμa*,δμs*)L22).
K((μ̅a,μ̅s)+ω(δμa*,δμs*))=K(μ̅a,μ̅s)ωR(μ̅a,μ̅s)L22+O(δμa*,δμs*L22).
K((μ̅a,μ̅s)+ω(δμa*,δμs*))=K(μ̅a,μ̅s)ωβR'̃(μ̅a,μ̅s)R(μ̅a,μ̅s)L22
+ O(δμa*,δμs*L22),
Sn1Ĩ(ŝ)[I(ŝ)]dŝ=Sn1[Ĩ(ŝ)]I(ŝ)dŝ,
Ĩ,gFPI=Sn1Ĩ(ŝ)[a1Δŝ(a2Δŝ)1I(ŝ)]dŝ.
Ĩ,gFPI=Sn1[(a2Δŝ)J̃(ŝ)][a1ΔŝJ(ŝ)]dŝ
=Sn1[a1ΔŝJ̃(ŝ)][(a2Δŝ)J(ŝ)]dŝ
=Sn1[a1Δŝ(a2Δŝ)1Ĩ(ŝ)]I(ŝ)dŝ
=gFPI,Ĩ.

Metrics