OSA's Digital Library

Optics Express

Optics Express

  • Editor: J. H. Eberly
  • Vol. 7, Iss. 13 — Dec. 18, 2000
  • pp: 492–506
« Show journal navigation

Scattering and absorption transport sensitivity functions for optical tomography

Oliver Dorn  »View Author Affiliations


Optics Express, Vol. 7, Issue 13, pp. 492-506 (2000)
http://dx.doi.org/10.1364/OE.7.000492


View Full Text Article

Acrobat PDF (303 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Optical tomography is modelled as an inverse problem for the time-dependent linear transport equation. We decompose the linearized residual operator of the problem into absorption and scattering transport sensitivity functions. We show that the adjoint linearized residual operator has a similar physical meaning in optical tomography as the ‘backprojection’ operator in x-ray tomography. In this interpretation, the geometric patterns onto which the residuals are backprojected are given by the same absorption and scattering transport sensitivity functions which decompose the forward residual operator. Moreover, the ‘backtransport’ procedure, which has been introduced in an earlier paper by the author, can then be interpreted as an efficient scheme for ‘backprojecting’ all (filtered) residuals corresponding to one source position simultaneously into the parameter space by just solving one adjoint transport problem. Numerical examples of absorption and scattering transport sensitivity functions for various situations (including applications with voids) are presented.

© Optical Society of America

1 The transport equation in optical tomography

We model optical tomography by the time-dependent linear transport equation

1cut+θ·u(x,θ,t)+(a(x)+b(x))u(x,θ,t)b(x)Sn1η(θ·θ)u(x,θ,t)dθ
=δ(t)δ(xxs)δ(θθs)inΩ×Sn1×[0,T].
(1)

In order to uniquely determine a solution of (1), we augment it by a homogeneous initial condition,

u(x,θ,0)=0inΩ×Sn1,
(2)

and an absorbing boundary condition,

u(x,θ,t)=0onΓ.
(3)

Here, Ω is a finite domain in ℝn(n=2,3) with smooth boundary ∂Ω. S n-1 denotes the set of direction vectors θ which have unit length in ℝn. v(x) denotes the outward unit normal to Ω at the point xΩ. The solution of (1)–(3), u(x, θ,t), describes the density of particles (photons) which travel in Ω at time t through the point x in the direction θ. The velocity c of the particles is assumed to be normalized to c=1cms-1, such that we will drop it in the notation from now on.

The symbol Γ- (resp. Γ+) denotes the set of points (x,θ,t) which correspond to incoming (resp. outgoing) radiation (particles) at the boundary Ω:

Γ±:={(x, θ,t)∊Ω×Sn -1×[0,T], ± ν(xθ>0}.

a(x) is the absorption cross section, b(x) the scattering cross section, and µ(x)=a(x)+b(x) the total cross section or attenuation. These parameters are assumed to be real and strictly positive, and they depend only on the position x. The absorption mean free path la is defined by la(x):=a -1(x), the scattering mean free path ls by ls(x):=b -1(x), and the transport mean free path l by l(x):=µ -1(x). Typical values in optical tomography are la≈1.0–10.0 cm, lbl≈0.005-0.01 cm [3

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

, 13

13. J. P. Kaltenbach and M. Kaschke, “Frequency and Time-Domain Modelling of Light Transport in Random Media,” in:G. Müller (ed.): Medical Optical Tomography: Functional Imaging and Monitoring (1993)

].

The scattering function η(θ·θ′) describes the probability for a particle which enters a scattering event with the direction of propagation θ′, to leave this event with the direction θ. It is normalized to

Sn1η(θ·θ)=1,,

which expresses particle conservation in pure scattering events. The dot-product in the argument indicates that η depends only on the cosine of the scattering angle θ·θ′=:cosϑ and, in particular, is independent of the location of the scattering event x. A frequently used scattering function is the following Henyey-Greenstein scattering function

η(θ·θ)=1g22(1+g22gcosϑ)3/2
(4)

with cos ϑ=θ·θ′, where g∊]-1,1[is some parameter. A value of g=0 means isotropic scattering, whereas values between g=0.9 and g=0.95 are more typical for situations in optical tomography and describe scattering events which are primarily forward directed. We will use this scattering function in our numerical experiments. The right hand side of (1) models the source distribution

q(x,θ,t)=δ(t)δ(xxs)δ(θθs).
(5)

It is a delta-like laser pulse which is transmitted at time t=0 at the position xs into the direction θs. The symbols δ(x-xs), δ(t) and δ(θ-θs) denote Dirac delta functions in the corresponding variables. If the source is positioned at the boundary, xsΩ, we assume that v(xsθs<0. In our numerical experiments in section 5 we will use v(xsθs=-1 (normally incident radiation).

The initial condition (2) indicates that there are no photons travelling in Ω at the starting time of our experiment. The boundary condition (3) indicates that during the experiment no photons enter the domain Ω from the outside. All photons inside of Ω originate from the source q which, however, might be situated at the boundary Ω. We want to mention that an alternative way of describing laser sources at the boundary is to use an inhomogeneous boundary condition instead of (3) and a zero source q=0. Both choices are equivalent, although care has to be taken when transforming one into the other. For more details see for example [5

5. K. M. Case and P. F. Zweifel, “Linear Transport Theory” (Plenum Press, New York, 1967)

].

A typical form of the measurement functional in optical tomography is given by

Ma,b(xr,tr)=S+n-1ν(xr)·θu(xr,θ,tr)onΩ×[0,T],
(6)

where u is the solution of (1)–(3) with source q and parameters a, b and xr and tr denote the receiver location and receiving time, respectively. The symbol Sn -1 + denotes the subset of direction vectors θSn -1 for which v(xθ>0. The inverse transport problem in optical tomography can now be formulated as follows.

Inverse transport problem of optical tomography. Assume that we measure for p≥1 different sources of the form (5) the corresponding data which are given by (6), where u solves the transport problem (1)–(3). Provided with this information, and knowing the scattering function η, we want to reconstruct both coefficient functions a(x) and b(x) inside of Ω simultaneously.

2 The linearized residual operator and its adjoint

In this and the following section we will assume that the source (5) is given and fixed. In section 4 we will indicate possible ways to extend these results to the more general case when data corresponding to more than one source positions are given.

We denote the space of parameters a or b by P, and the space of data Ma,b by D. The residuals R(a, b) are defined as the difference between the physically measured data (x,t) and the data which correspond to the parameter distribution (a, b). More formally we have with (6)

R:P×PD,R(a,b)(xr,tr)=Ma,b(xr,tr)G˜(xr,tr).
(7)

The operator (7) is nonlinear in (a, b), since the solution u in (6) depends in a nonlinear way on (a, b). When solving the inverse problem, we want to find a set of parameters (â,) from the (usually noisy) data such that the residuals (7) are minimized in some given norm.

Since the operator R is nonlinear, many reconstruction approaches involve the calculation of its ‘derivative’ or ‘linearized operator’ R′a b which has to be calculated at the most recent best guess (a, b) for the parameters. This operator is often called ‘Fréchet-derivative’, ‘Jacobian’, or ‘sensitivity matrix’. It is a mapping from the parameter space P×P into the data space D.

A possible way of deriving explicit forms for the linearized residual operator is to perturb the parameter functions a(x)→a(x)+δ a(x) and b(x)→b(x)+δb(x) and plug this into (1)–(3). Assuming that the corresponding solution u(x) of (1)–(3) responds to this perturbation according to u(x)→u(x)+w(x), and neglecting all terms in (1) which are of higher than linear order in δa, δb and w, we arrive at the following result [9

9. O. Dorn, “Das inverse Transportproblem in der Lasertomographie,” Thesis, Preprints “Angewandte Mathematik und Informatik” 7/97-N, Münster (1997)

, 19

19. F. Natterer, “Numerical Solution of Bilinear Inverse Problems,” Preprints “Angewandte Mathematik und Informatik” 19/96-N, Münster (1996)

].

Expression for the linearized residual operator. The linearized residual operator R′a,b is given by

Ra,b:P×PD,Ra,b(δa,δb)(xr,tr)=S+n1ν(xr)·θω(xr,θ,tr)
(8)

where w solves the linearized equation

ωt+θ·ω(x,θ,t)+(a(x)+b(x))ω(x,θ,t)b(x)Sn1η(θ·θ)ω(x,θ,t))
=Qδa(x,θ,t)+Qδb(x,θ,t)
(9)

in Ω×Sn -1×[0,T] with the homogeneous initial condition

ω(x,θ,0)=0inΩ×Sn1,
(10)

and an absorbing boundary condition

ω(x,θ,t)=0onΓ.
(11)

Here, thescattering sourcesQδa and Qδb are defined as

Qδa(x,θ,t)=δa(x)u(x,θ,t)
(12)

and

Qδb(x,θ,t)=δb(x)(u(x,θ,t)Sn1η(θ·θ)u(x,θ,t))
(13)

where u is the solution of (1)–(3) with parameters a and b.

Notice that, for a given function pair (δa(xsc), δb(xsc)), where the index sc stands for ‘scatter’, the value of R′a,b(δa, δb) is a function in the variables xr and tr, where xr is the receiver location and tr the receiver time. This explains the somewhat complicated notation in (8). The physical interpretation of this result is that the perturbations δa and δb create scattering sources Qδa and Qδb, which give rise to a distribution w(x, θ, t) (which can be positive or negative) of virtual ‘secondary particles’ propagating in the unperturbed medium to the receivers, where they are detected as the (linearized) residuals in the data.

The adjoint linearized residual operator R′*a,b is formally defined by the identity

Ra,b(δa,δb),ςD=(δa,δb),Ra,b*ςP×P
(14)

where ζ is some vector in the data space D, and where the brackets denote the inner products in the underlying function spaces P×P and D. The adjoint linearized residual operator is a mapping from the data space D into the parameter space P×P.

An explicit expression for the action of the adjoint linearized residual operator on an arbitrary vector ζ of the data space can be derived by using Green’s formula for the linear transport equation in time domain. Next, we formulate this result, and refer for a derivation of it to [9

9. O. Dorn, “Das inverse Transportproblem in der Lasertomographie,” Thesis, Preprints “Angewandte Mathematik und Informatik” 7/97-N, Münster (1997)

, 10

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

].

Expression for the adjoint linearized residual operator. Let z denote the solution of the following adjoint linear transport equation

ztθ·z(x,θ,t)+(a(x)+b(x))z(x,θ,t)b(x)Sn1η(θ·θ)z(x,θ,t)
=0inΩ×Sn1×[0,T]
(15)

with thefinal value condition

z(x,θ,T)=0inΩ×Sn1,
(16)

and the inhomogeneous outgoing boundary condition

z(x,θ,t)=ς(x,t)uniformlyinθonΓ+,
(17)

and let u be the solution of the forward problem (1)–(3). Then we have

(Ra,b*ς)(x)=(Δa(x),Δb(x))
(18)

with

Δa(x)=[0,T]Sn1u(x,θ,t)z(x,θ,t)dt,
(19)
Δa(x)=[0,T]Sn1[Sn1η(θ·θ)u(x,θ,t)u(x,θ,t)]z(x,θ,t)dt.
(20)

Solving (15)–(17) is calledbacktransport’.

We mention that also the adjoint linearized residual operator admits a nice physical interpretation. The argument ζ of this operator is an element of the data space D, i.e., it is a function of detector position xr and detection time tr. The values of ζ, which in our applications will be the (filtered) residuals, are attached as time-dependent adjoint sources at these receiver positions in (17), and transported backward in direction and in time into the medium Ω by solving (15)–(17). Notice the sign change in front of the time derivative t and the space derivative θ·∇ in (15) compared to (1), which indicates a reversed time and direction. This is also the reason why we have used a final value condition (16) and an ‘outgoing flux condition’ (17) for uniquely specifying the solution of (15). Notice also that the (filtered) residuals ζ are applied uniformly in all directions in (15), which compensates for the averaging θ-integration of the measurement process (6).

Loosely speaking, these virtual ‘adjoint particles’ are trying to trace back the paths of those virtual particles which have been created by the parameter perturbation (δa, δb) and which have caused the mismatch in the data. The correction to our best guess (which is just R′*a,bζ as it is described in section 4) is then calculated by combining these backpropagated densities with the actual densities u of our forward problem (1)–(3). For example, if no particles of the forward solution reach a given point xsc during the experiment, i.e. u(xsc, θ, t)=0 for all θSn -1 and all t∊[0,T], then the update (R′*a,bζ)(x) in (19), (20) will be zero at this location xsc. This makes sense physically since no secondary particles could have been generated in (12), (13) at this point in the medium by a parameter change (δa(xsc), δb(xsc)) and reach the detector via (9). This means also that the ‘sensitivity’ of our source-receiver pair to parameter changes at this location will be zero. This observation motivates the introduction of (linearized) ‘sensitivity functions’, which quantify the sensitivity of a given source-receiver pair to parameter perturbations at each point xsc in the medium. We will discuss these sensitivity functions for the linear transport equation in the following section.

3 Transport sensitivity functions

Transport sensitivity functions. For a given source (5), we can find functions ψa(xr,tr;xsc) and ψb(xr,tr;xsc) with the following properties.

1.)‘ProjectionStep: The action of the linearized residual operator R′a,b on the perturbation ((δa, δb) in parameter space can be described as follows

(Ra,b(δa,δb)(xr,tr)=
ΩdxscΨa(xr,tr;xsc)δa(xsc)+ΩdxscΨb(xr,tr;xsc)δb(xsc)
(21)

2.) ‘BackprojectionStep: The action of the adjoint linearized residual operator R′*a,b on the vector ζ in data space can be described as follows

(Ra,b*ς)(xsc)=
(Ωdxr[0,T]dtrΨa(xr,tr;xsc)ς(xr,tr),Ωdxr[0,T]dtrΨb(xr,tr;xsc)ς(xr,tr))
(22)

The functions ψa(xr,tr;xsc) and ψb(xr,tr;xsc) are calledtransport sensitivity functions’.

Formula (21) says that a parameter perturbation (δa(xsc), δb(xsc)) at a position xsc with a large sensitivity value ψa,b(xr,tr;xsc) will have a relatively strong influence on the data measured at the location xr and at the time tr. Formula (22), on the other hand, says that a residual at a receiver with position xr, which is detected at time tr, will produce relatively large updates of the parameters (a, b) in the medium at positions xsc where the corresponding sensitivity functions ψa,b(xr,tr;xsc) have large values.

In our case, we want to solve a nonlinear inverse problem. The adjoint linearized residual operator R′*a,b changes here (as part of an iterative reconstruction scheme) with the latest best guess (a, b). Moreover, it is not applied to the data but to the differences between the data and the calculated data. Moreover, the geometric patterns onto which the data are ‘backprojected’ are different here. Instead of backprojecting over lines, as it is done in x-ray tomography, we have to backproject onto more complicated shapes which are given by the sensitivity functions ψa,b(xr,tr;xsc). These sensitivity functions depend on the source and receiver locations, the most recent parameter distribution in the medium, and the detection time of the given data.

4 The Transport-BackTransport (TBT) algorithm

The ‘backtransport’ procedure described in section 2 provides us with a very efficient scheme to ‘backproject’ all given data (residuals) ζ corresponding to one source position, and for the most recent parameter distribution (a, b), simultaneously by just solving one adjoint transport problem (15)–(17) on a computer. In many applications, however, we have data given which correspond to many source positions. A possible way of addressing these problems and making use of the adjoint scheme is to consider only the data for one source at a time. We want to outline this procedure in the following.

When employing a linearization (Newton-type) approach to a data set corresponding to one given source position, we can calculate a correction (δa, δb) to our best guess (a, b) by looking for a solution of the linearized problem

Ra,b(δa,δb)=R(a,b).
(23)

Because of the limited amount of data which are used in this step, we assume that (23) is underdetermined. The (generalized) solution of this problem has the well-known form

(δa,δb)=Ra,b*(Ra,bRa,b*)1R(a,b).
(24)

Let us assume that we can find a good approximation C to the unknown operator (R′a,b R′*a,b)-1, which maps from the data space to the data space and which can be considered as a ‘filtering operator’. We have shown that all we have to do in order to find (δa, δb) is to solve one forward transport problem (1)–(3) with our given source q, and one adjoint transport problem (15)–(17) with the (filtered) residual values ζ=CR(a, b) attached as ‘adjoint sources’ to the receiver positions. The correction (δa, δb) is then calculated from these two solutions u and z by the formulas (19) and (20). We mention that additional regularization criteria can be incorporated in (24), but these refinements of the inversion scheme will not be discussed in the present paper.

Fig. 1. The geometries of the four experiments. Displayed is the distribution of the scattering cross section b. The values are b=1.0 cm-1 inside the voids, and b=100.0 cm-1 in the background. Top left: experiment 1; Top right: experiment 2; Bottom left: experiment 3; Bottom right: experiment 4. Indicated in the top left image are also the boundary points P1, P2 and P3, where the sources and receivers are located.

5 Numerical examples of transport sensitivity functions

In our numerical experiments we use a finite-differences discretization of the linear transport equation (1)(3) in 2D which has been described in [9

9. O. Dorn, “Das inverse Transportproblem in der Lasertomographie,” Thesis, Preprints “Angewandte Mathematik und Informatik” 7/97-N, Münster (1997)

, 10

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

]. Our computational domain Ω is a square of the size 6×6 cm2 which is discretized into 120×120 elementary cells or pixels. One individual time step lasts 0.2 s. Assuming a (normalized) velocity of c=1 cm/s, a photon needs for example 6 s or 30 time steps in order to travel from one boundary of the domain parallel to the axes to the opposite boundary without being scattered. The directional variable θ is discretized into 12 individual direction vectors equidistantly distributed over the unit circle S 1, four of them pointing in the directions of the positive and negative x and y coordinate axes. For a more detailed description of this discretization scheme for (1)–(3) we refer to [9

9. O. Dorn, “Das inverse Transportproblem in der Lasertomographie,” Thesis, Preprints “Angewandte Mathematik und Informatik” 7/97-N, Münster (1997)

, 10

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

].

We consider four different physical situations, which we call experiments 1–4 and which are visualized in Figure 1. The first experiment considers a homogeneous parameter distribution with b=100 cm-1 and a=0.1 cm-1 everywhere. In the other three experiments, we have embedded into this background distribution so-called ‘clear regions’ or ‘voids’ of various shapes. These clear regions can be found in many physical situations where near-infrared light is used for medical imaging [4

4. S. R. Arridge, H. Dehgani, M. Schweiger, and E. Okada, “The Finite Element Model for the Propagation of Light in Scattering Media: A Direct Method for Domains with Non-Scattering Regions,” Medical Physics , 27 (1), 252–264 (2000) [CrossRef] [PubMed]

, 15

15. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data” Inverse Problems 15, 1375–1391 (1999) [CrossRef]

, 20

20. E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36 (1), 21–31 (1997) [CrossRef] [PubMed]

, 21

21. J. Riley, H. Dehghani, M. Schweiger, S. R. Arridge, J. Ripoll, and M. Nieto-Vesperinas, “3D Optical Tomography in the Presence of Void Regions,” to appear in Optics Express, (this issue)

]. Therefore, it is of practical importance to understand the behavior of the forward and the adjoint transport problems in the presence of voids. In our three experiments including voids the parameters inside the clear regions are chosen to be b=1.0 cm-1 and a=0.01 cm-1. Their shapes can be seen in Figure 1. We mention that in all of our numerical examples a Henyey-Greenstein scattering law (4) is used with parameter g=0.9.

Fig. 2. Left: Data corresponding to a source at boundary point P1=(40, 0) and a receiver at P3=(60, 120). Right: Data corresponding to a source at boundary point P2=(120, 40) and a receiver at P3=(60, 120). In both images the curves correspond to (from bottom to top) experiment 1 (blue, dash-dotted), 4 (red, solid), 2 (green, stars), and 3 (magenta, dashed). The vertical lines in the figures indicate at which time steps (T1=50 and T2=120) the sensitivity functions displayed in this paper have been calculated.

In Figure 2 we have displayed some typical data Ma,b(xr, tr), defined in (6), for these four experiments. In the first image, a delta-like laser source (5) is positioned at the boundary point P1(which is marked in Figure 1 in the upper left image) with coordinates (40, 0). The receiver is situated at the boundary point P3 (see again Figure 1) and has the coordinates (60, 120). In the second image, the source is positioned at the boundary point P2 with coordinates (120,40), and the receiver is again situated at the boundary point P3. The data curves correspond in both images, from bottom to top, to the experiments 1 (blue, dash-dotted), 4 (red, solid), 2 (green, stars), and 3 (magenta, dashed). We see that the presence of clear regions has a significant influence on the first arrival times as well as on the amplitudes and shapes of the data.

In particular it can be noticed that the clear layers act as a ‘wave guide’ for the early particles. In addition, during the time of our experiments, much more particles actually reach the two detectors considered here than in the case without clear regions. Comparing the data of experiments 2 and 3 (with a source located at P2) on the right hand side of Figure 2, we see that for very early times the two curves coincide, but approximately at time step 40 they split. This is probably due to the disc-like void region in experiment 3 which allows particles to propagate more quickly than in the strongly scattering background of experiment 2. These effects would not be visible if a purely diffusion-based model would have been used instead of the transport model for modelling the light propagation in these two experiments. We mention, however, that a hybrid diffusion-radiosity based model has been reported recently in [4

4. S. R. Arridge, H. Dehgani, M. Schweiger, and E. Okada, “The Finite Element Model for the Propagation of Light in Scattering Media: A Direct Method for Domains with Non-Scattering Regions,” Medical Physics , 27 (1), 252–264 (2000) [CrossRef] [PubMed]

, 21

21. J. Riley, H. Dehghani, M. Schweiger, S. R. Arridge, J. Ripoll, and M. Nieto-Vesperinas, “3D Optical Tomography in the Presence of Void Regions,” to appear in Optics Express, (this issue)

] which should provide an interesting alternative to the pure transport model used here for the modelling of scattering media with voids.

In the Figures 39 we have displayed various typical absorption sensitivity functions ψa(xr, tr; xsc) and scattering sensitivity functions ψb(xr, tr; xsc) for these four experiments. (Notice that the figures for the absorption sensitivity functions actually show -ψa(xr, tr; xsc) for improving the presentation.) The images are generally arranged in the same way as in Figure 1. In order to simplify the comparison of these figures with each other, we have normalized these functions, using the MATLAB routine norm(.,‘fro’) for calculating the normalization factors (Frobenius norms) N. The corresponding normalization factors are given in the captions of the figures.

Fig. 3. Transport absorption sensitivity functions -ψa(xr, tr; xsc) for a source located at P1 and a receiver at P3. The receiving time is tr=10 s, which corresponds to the time step T1=50. Top left: experiment 1 (N=99.9); Top right: experiment 2 (N=664.1); Bottom left: experiment 3 (N=1000.0); Bottom right: experiment 4 (N=221.8).

Figure 3 displays -ψa(xr, tr; xsc) where the source is situated at position P1 and the receiver location xr is at P3. The receiving time is tr=10 s, which corresponds to the time step T1=50. Therefore, looking at the data curves on the left hand side of Figure 2, this figure describes the sensitivity of the early photons (which are sometimes called ‘snake’ photons) to parameter perturbations in the imaging domain. We can clearly see the ‘wave-guiding effect’ which occurs in the experiments 2 and 3 due to the near-boundary clear layer.

Figure 4 shows the same situation as Figure 3, but with a receiving time of tr=24 s, which corresponds to time step T2=120. See again the image on the left hand side of Figure 2 for locating these receiver times in the general form of the data curves. We see that the shapes of the sensitivity functions at these later times are much broader than those of the snake photons. In addition, the wave-guiding effect of the clear boundary layers has disappeared. Figures 5 and 6 show the sensitivity functions which correspond to Figures 3 and 4 but where the source is now located at position P2 instead of P1. The discussion of these two figures is completely analogous to the discussion of Figures 3 and 4 and is therefore omitted here.

Figures 7 and 9 show the scattering sensitivity functions ψb(xr, tr; xsc) for the same situations as in Figures 3 and 4. We see that the data corresponding to early times as well as to later times are equally sensitive to parameter changes inside the clear layers. Notice also that the norms N of these sensitivity functions are much smaller than those for the absorption parameter. This is due to the strong scattering in the domain. Directional information corresponding to local perturbations in the scattering parameter is quickly lost due to the many succeeding scattering events of the photons.

Fig. 4. Transport absorption sensitivity functions -ψa(xr, tr; xsc) for a source located at P1 and a receiver at P3. The receiving time is tr=24 s, which corresponds to the time step T2=120. Top left: experiment 1 (N=3.0×104); Top right: experiment 2 (N=5.4×104); Bottom left: experiment 3 (N=7.3×104); Bottom right: experiment 4 (N=4.1×104).
Fig. 5. Transport absorption sensitivity functions -ψa(xr, tr; xsc) for a source located at P2 and a receiver at P3. The receiving time is tr=10 s, which corresponds to the time step T1=50. Top left: experiment 1 (N=5.6×103); Top right: experiment 2 (N=8.5×104); Bottom left: experiment 3 (N=1.0×105); Bottom right: experiment 4 (N=1.3×104).
Fig. 6. Transport absorption sensitivity functions -ψa(xr, tr; xsc) for a source located at P2 and a receiver at P3. The receiving time is tr=24 s, which corresponds to the time step T2=120. Top left: experiment 1(N=9.8×104); Top right: experiment 2(N=1.7×105); Bottom left: experiment 3(N=2.3×105); Bottom right: experiment 4 (N=1.4×105).

In Figure 8 we have displayed ψb(xr, tr; xsc) as in example 3 of Figure 7, but calculated with different numerical discretizations. Consult the caption of Figure 8 for more details. We see that the structures inside the void regions are reproduced reliably (with small numerical fluctuations), which indicates that these structures are most likely not caused by numerical noise but that they represent a meaningful physical phenomenon. It would be desirable, however, to compare these results with a completely independent numerical model in order to confirm these observations.

Finally, we want to mention that the theory presented in this paper can be generalized to situations where additional physical properties of the probing light is taken into account. One example would be to use polarized light [6

6. S. Chandrasekhar, “Radiative Transfer” (Dover, New York, 1960)

, 23

23. J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31, 6535–6546 (1992) [CrossRef] [PubMed]

], which leads to a transport equation for a 2×2 coherence matrix describing light intensity and polarization. The numerical modelling of the propagation of polarized light in scattering media is, however, more difficult and is a topic of ongoing research [16

16. M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissues,” to appear in J. Opt. Soc. Am. (2000)

, 22

22. L. Ryzhik, G. Papanicolaou, and J. B. Keller, “Transport equations for elastic and other waves in random media,” Wave Motion 24, 327–370 (1996) [CrossRef]

].

6 Appendix

In this section we give explicit expressions for the transport sensitivity functions, and indicate how to calculate them practically. We define the transport Green functions G[xsc, θsc,tsc](x, θ, t) as solutions of the problem

Gt+θ·G(x,θ,t)+(a(x)+b(x))G(x,θ,t)b(x)Sn1η(θ·θ)G(x,θ,t)
=δ(xxsc)δ(θθsc)δ(ttsc)
(25)

in Ω×Sn -1×[0,T] with the homogeneous initial condition

G(x,θ,0)=0inΩ×sn1,
(26)

and the absorbing boundary condition

G(x,θ,t)=0onΓ.
(27)
Fig. 7. Transport scattering sensitivity functions ψb(xr, tr; xsc) for a source located at P1 and a receiver at P3. The receiving time is tr=10 s, which corresponds to the time step T1=50. Top left: experiment 1 (N=1.3); Top right: experiment 2 (N=166.6); Bottom left: experiment 3 (N=164.3); Bottom right: experiment 4 (N=3.45).
Fig. 8. Displayed is the same transport scattering sensitivity function as in example 3 of Figure 7, but calculated with different numerical discretizations. Left: 12 direction vectors are used as in Figure 7, but a finer spatial grid of 240×240 pixels instead of 120×120 pixels. Right: A 120×120 grid is used as in Figure 7, but 24 discretized directions instead of 12 directions. Shown are ψb(xr, tr; xsc) for a source located at P1 and a receiver at P3. The receiving time step is T1=50.
Fig. 9. Transport scattering sensitivity functions ψb(xr, tr; xsc) for a source located at P1 and a receiver at P3. The receiving time is tr=24 s, which corresponds to the time step T2=120. Top left: experiment 1 (N=144.3); Top right: experiment 2 (N=443.0); Bottom left: experiment 3 (N=619.1); Bottom right: experiment 4 (N=204.0).

The index sc indicates an arbitrary point in the scattering region. The adjoint transport Green functions Ĝ[xr, θr, tr](x, θ,t) are defined accordingly by

Ĝtθ·Ĝ(x,θ,t)+(a(x)+b(x))Ĝ(x,θ,t)b(x)Sn1η(θ·θ)Ĝ(x,θ,t)
=δ(xxr)δ(θθr)δ(ttr)
(28)

in Ω×Sn -1×[0,T] with the ‘final value condition’

Ĝ(x,θ,T)=0inΩ×Sn1,
(29)

and the outgoing boundary condition

Ĝ(x,θ,t)=0onΓ+.
(30)

The functions Φa[xr, θr, tr] are defined by

Φa[xr,θr,tr](xsc,θsc,tsc)=Ĝ[xr,θr,tr](xsc,θsc,tsc)G[xs,θs,0](xsc,θsc,tsc).
(31)

Let us assume for the moment that we only have a perturbation in the absorption parameter, which means δb=0. Then the linearized residual operator corresponding to our source (5) can be written in the explicit form

(Ra,b(δa,0))(xr,tr)=-S+n1rΩdxscSn1sc[0,T]dtscν(xr)·θr
×G[xs,θs,0](xsc,θsc,tsc)G[xsc,θsc,tsc](xr,θr,tr)δa(xsc)
(32)

where the index s stands for ‘source’, and r stands for ‘receiver’. Using the following generalized reciprocity relation

G[xsc,θsc,tsc](xr,θr,tr)=Ĝ[xr,θr,tr](xsc,θsc,tsc),
(33)

and defining the absorption sensitivity functions ψa(xr,tr;xsc) by

Ψa(xr,tr;xsc)=S+n1rν(xr)·θr[0,T]dtscSn1scΦa[xr,θr,tr](xsc,θsc,tsc),
(34)

we can write (32) in the shorter form

(Ra,b(δa,0))(xr,tr)=ΩdxscΨa(xr,tr;xsc)δa(xsc),
(35)

which is the first term in (21).

In order to derive the second term in (21), we consider now a perturbation (0, δb) ∊ P×P and define the functions Φb[xr, θr, tr] by

Φb[xr,θr,tr](xsc,θsc,tsc)=Ĝ[xr,θr,tr](xsc,θsc,tsc)×
[G[xs,θs,0](xsc,θsc,tsc)Sn1η(θsc·θ)G[xs,θs,0](xsc,θ,tsc)].
(36)

Introducing the scattering sensitivity functions ψb(xr,tr;xsc) as

Ψb(xr,tr,xsc)=S+n1rν(xr).θr[0,T]dtscSn1scΦb[xr,θr,tr](xsc,θsc,tsc),
(37)

we arrive with the same arguments as before at the second term in (21).

The corresponding expressions for the adjoint linearized residual operator R′*a,b are derived in a similar way, such that we will not present them here explicitly. In the derivation we make use of the fact that a boundary condition z(x, θ, t)=ζ(x, t) on Γ+ can be replaced by an equivalent ‘surface source distribution’ q∂Ω=v(x)·θζ(x, t) along the boundary. For more details see for example [5

5. K. M. Case and P. F. Zweifel, “Linear Transport Theory” (Plenum Press, New York, 1967)

].

Practically, we only have to solve one forward transport problem (1)–(3) and one adjoint transport problem (15)–(17) with ζ,(x,t)=δ(x-xr)δ(t-tr) in order to calculate the sensitivity functions ψa(xr, tr; xsc) and ψb(xr, tr; xsc) corresponding to a given receiver position xr and a receiver time tr. The sensitivity functions are then automatically given by ψa(xr,tr; xsc)=Δa(xsc) and ψb(xr, tr; xsc)=Δb(xsc) with Δa and Δb defined in (19) and (20).

7 Acknowledgment

This work was supported by the NSF grant DMS97-09320 and the NSERC CRD Grant 80357. The article benefitted from helpful comments by Miguel Moscoso.

References and links

1.

S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. 34 (34), 7395–7409 (1995) [CrossRef] [PubMed]

2.

S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part II: Finite-element-method calculations,” Appl. Opt. 34 (31), 8026–8037 (1995) [CrossRef] [PubMed]

3.

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

4.

S. R. Arridge, H. Dehgani, M. Schweiger, and E. Okada, “The Finite Element Model for the Propagation of Light in Scattering Media: A Direct Method for Domains with Non-Scattering Regions,” Medical Physics , 27 (1), 252–264 (2000) [CrossRef] [PubMed]

5.

K. M. Case and P. F. Zweifel, “Linear Transport Theory” (Plenum Press, New York, 1967)

6.

S. Chandrasekhar, “Radiative Transfer” (Dover, New York, 1960)

7.

S. B. Colak, D. G. Papanioannou, G. W. ’t. Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens, J. B. M. Melissen, and N. A. A. J. van Asten, “Tomographic image reconstruction from optical projections in light-diffusing media,” Appl. Opt. 36, 180–213 (1977) [CrossRef]

8.

T. Dierkes, “Rekonstruktionsverfahren zur optischen Tomographie,” Thesis, Preprints “Angewandte Mathematik und Informatik” 16/00-N, Münster (2000)

9.

O. Dorn, “Das inverse Transportproblem in der Lasertomographie,” Thesis, Preprints “Angewandte Mathematik und Informatik” 7/97-N, Münster (1997)

10.

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

11.

S. Feng, Z. A. Zeng, and B. Chance, “Photon migration in the presence of a single defect: a perturbation analysis,” Appl. Opt. 34, 3826–3837 (1995) [CrossRef] [PubMed]

12.

R. J. Gaudette, C. A. Brooks, M. E. Kilmer, E. L. Miller, T. Gaudette, and D. Boas, “A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient,” Phys. Med. Biol. 45 (4) 1051–1070 (2000) [CrossRef] [PubMed]

13.

J. P. Kaltenbach and M. Kaschke, “Frequency and Time-Domain Modelling of Light Transport in Random Media,” in:G. Müller (ed.): Medical Optical Tomography: Functional Imaging and Monitoring (1993)

14.

M. V. Klibanov, T. R. Lucas, and R. M. Frank, “A fast and accurate imaging algorithm in optical/diffusion tomography,” Inverse Problems 13, 1341–1361 (1997) [CrossRef]

15.

V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data” Inverse Problems 15, 1375–1391 (1999) [CrossRef]

16.

M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissues,” to appear in J. Opt. Soc. Am. (2000)

17.

F. Natterer, “The Mathematics of Computerized Tomography” (Stuttgart, Teubner, 1986)

18.

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

19.

F. Natterer, “Numerical Solution of Bilinear Inverse Problems,” Preprints “Angewandte Mathematik und Informatik” 19/96-N, Münster (1996)

20.

E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36 (1), 21–31 (1997) [CrossRef] [PubMed]

21.

J. Riley, H. Dehghani, M. Schweiger, S. R. Arridge, J. Ripoll, and M. Nieto-Vesperinas, “3D Optical Tomography in the Presence of Void Regions,” to appear in Optics Express, (this issue)

22.

L. Ryzhik, G. Papanicolaou, and J. B. Keller, “Transport equations for elastic and other waves in random media,” Wave Motion 24, 327–370 (1996) [CrossRef]

23.

J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31, 6535–6546 (1992) [CrossRef] [PubMed]

24.

J. C. Schotland, J. C. Haselgrove, and J. S. Leigh, “Photon hitting density,” Appl. Opt. 32, 448–453 (1993) [CrossRef] [PubMed]

OCIS Codes
(170.3880) Medical optics and biotechnology : Medical and biological imaging
(170.5280) Medical optics and biotechnology : Photon migration
(170.6960) Medical optics and biotechnology : Tomography
(170.7050) Medical optics and biotechnology : Turbid media

ToC Category:
Focus Issue: Diffuse optical tomography

History
Original Manuscript: October 27, 2000
Published: December 18, 2000

Citation
Oliver Dorn, "Scattering and absorption transport sensitivity functions for optical tomography," Opt. Express 7, 492-506 (2000)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-7-13-492


Sort:  Journal  |  Reset  

References

  1. S. R. Arridge, "Photon-measurement density functions. Part I: Analytical forms," Appl. Opt. 34 (34), 7395-7409 (1995) [CrossRef] [PubMed]
  2. S. R. Arridge and M. Schweiger, "Photon-measurement density functions. Part II: Finite-element-method calculations," Appl. Opt. 34 (31), 8026-8037 (1995) [CrossRef] [PubMed]
  3. S. R. Arridge, "Optical tomography in medical imaging," Inverse Problems 15 (2), R41-R93 (1999) [CrossRef]
  4. S. R. Arridge, H. Dehgani, M. Schweiger and E. Okada, "The Finite Element Model for the Propagation of Light in Scattering Media: A Direct Method for Domains with Non-Scattering Regions," Medical Physics, 27 (1), 252-264 (2000) [CrossRef] [PubMed]
  5. K. M. Case and P. F. Zweifel, "Linear Transport Theory" (Plenum Press, New York, 1967)
  6. S. Chandrasekhar, "Radiative Transfer" (Dover, New York, 1960)
  7. S. B. Colak, D. G. Papanioannou, G. W. 't. Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens, J. B. M. Melissen and N. A. A. J. van Asten, "Tomographic image reconstruction from optical projections in light-diffusing media," Appl. Opt. 36, 180-213 (1977) [CrossRef]
  8. T. Dierkes, "Rekonstruktionsverfahren zur optischen Tomographie," Thesis, Preprints "Angewandte Mathematik und Informatik" 16/00-N, Münster (2000)
  9. O. Dorn, "Das inverse Transportproblem in der Lasertomographie," Thesis, Preprints "Angewandte Mathematik und Informatik" 7/97-N, Münster (1997)
  10. O. Dorn, "A transport-backtransport method for optical tomography," Inverse Problems 14, 1107-1130 (1998) [CrossRef]
  11. S. Feng, Z. A. Zeng and B. Chance, "Photon migration in the presence of a single defect: a perturbation analysis," Appl. Opt. 34, 3826-3837 (1995) [CrossRef] [PubMed]
  12. R. J. Gaudette, C. A. Brooks, M. E. Kilmer, E. L. Miller, T. Gaudette, D. Boas, "A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient," Phys. Med. Biol. 45 (4) 1051-1070 (2000) [CrossRef] [PubMed]
  13. J. P. Kaltenbach and M. Kaschke, "Frequency and Time-Domain Modelling of Light Transport in Random Media," in: Müller, G. (ed.): Medical Optical Tomography: Functional Imaging and Monitoring (1993)
  14. M. V. Klibanov, T. R. Lucas and R. M. Frank, "A fast and accurate imaging algorithm in optical/diffusion tomography," Inverse Problems 13, 1341-1361 (1997) [CrossRef]
  15. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen and J. P. Kaipio, "Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data" Inverse Problems 15, 1375-1391 (1999) [CrossRef]
  16. M. Moscoso, J. B. Keller and G. Papanicolaou, "Depolarization and blurring of optical images by biological tissues," to appear in J. Opt. Soc. Am. (2000)
  17. F. Natterer, "The Mathematics of Computerized Tomography" (Stuttgart, Teubner, 1986)
  18. F. Natterer and F. Wübbeling, "A propagation-backpropagation method for ultrasound tomography," Inverse Problems 11, 1225-1232 (1995) [CrossRef]
  19. F. Natterer, "Numerical Solution of Bilinear Inverse Problems," Preprints "Angewandte Mathematik und Informatik" 19/96-N, Münster (1996)
  20. E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope and D. T. Delpy, "Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head," Appl. Opt. 36 (1), 21-31 (1997) [CrossRef] [PubMed]
  21. J. Riley, H. Dehghani, M. Schweiger, S. R. Arridge, J. Ripoll and M. Nieto-Vesperinas, "3D Optical Tomography in the Presence of Void Regions," Opt. Express 7, 462-467 (2000), http://www.opticsexpress.org/oearchive/source/26894.htm
  22. L. Ryzhik, G. Papanicolaou and J. B. Keller, "Transport equations for elastic and other waves in random media," Wave Motion 24, 327-370 (1996) [CrossRef]
  23. J. M. Schmitt, A. H. Gandjbakhche and R. F. Bonner, "Use of polarized light to discriminate short-path photons in a multiply scattering medium," Appl. Opt. 31, 6535-6546 (1992) [CrossRef] [PubMed]
  24. J. C. Schotland, J. C. Haselgrove and J. S. Leigh, "Photon hitting density," Appl. Opt. 32, 448-453 (1993) [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited