## 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.

©2000 Optical Society of America

## 1 The transport equation in optical tomography

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

$$=\delta \left(t\right)\delta (x-{x}_{\mathit{s}})\mathit{\delta}(\mathit{\theta}-{\mathit{\theta}}_{\mathit{s}})\phantom{\rule{.9em}{0ex}}\mathrm{in}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\Omega \times {S}^{\mathit{n}-1}\times \left[0,T\right].$$

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

and an absorbing boundary condition,

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*)∊*∂*Ω×*S ^{n}*

^{-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 *l _{a}* is defined by

*l*(

_{a}*x*):=

*a*

^{-1}(

*x*), the scattering mean free path

*l*by

_{s}*l*(

_{s}*x*):=

*b*

^{-1}(

*x*), and the transport mean free path

*l*by

*l*(

*x*):=

*µ*

^{-1}(

*x*). Typical values in optical tomography are

*l*≈1.0–10.0 cm,

_{a}*l*≈

_{b}*l*≈0.005-0.01 cm [3, 13].

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

${\int}_{{S}^{n-1}}\eta (\theta \xb7{\theta}^{\prime})\mathit{d\theta}=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

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

It is a delta-like laser pulse which is transmitted at time *t*=0 at the position *x _{s}* into the direction

*θ*. The symbols

_{s}*δ*(

*x*-

*x*),

_{s}*δ*(

*t*) and

*δ*(

*θ*-

*θ*) denote Dirac delta functions in the corresponding variables. If the source is positioned at the boundary,

_{s}*x*∊

_{s}*∂*Ω, we assume that

*v*(

*x*)·

_{s}*θ*<0. In our numerical experiments in section 5 we will use

_{s}*v*(

*x*)·

_{s}*θ*=-1 (normally incident radiation).

_{s}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].

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

where *u* is the solution of (1)–(3) with source *q* and parameters *a*, *b* and *x _{r}* and

*t*denote the receiver location and receiving time, respectively. The symbol

_{r}*S*

^{n}^{-1}

_{+}denotes the subset of direction vectors

*θ*∊

*S*

^{n}^{-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 *M _{a,b}* by

*D*. The residuals

*R*(

*a*,

*b*) are defined as the difference between the physically measured data

*G̃*(

*x*,

*t*) and the data which correspond to the parameter distribution (

*a*,

*b*). More formally we have with (6)

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 (*â*,* b̂*) from the (usually noisy) data *G̃* 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, 19].

Expression for the linearized residual operator. *The linearized residual operator*
*R′ _{a,b}*

*is given by*

*where*
* w *
*solves the linearized equation*

$$={Q}_{\mathit{\delta a}}(x,\theta ,t)+{Q}_{\mathit{\delta b}}(x,\theta ,t)$$

in Ω×*S ^{n}*

^{-1}×[0,

*T*]

*with the homogeneous initial condition*

*and an absorbing boundary condition*

*Here, the* ‘*scattering sources*’ *Q _{δa}*

*and*

*Q*

_{δb}*are defined as*

and

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

Notice that, for a given function pair (*δa*(*x _{sc}*),

*δb*(

*x*)), where the index

_{sc}*sc*stands for ‘scatter’, the value of

*R′*(

_{a,b}*δa*,

*δb*) is a function in the variables

*x*and

_{r}*t*, where

_{r}*x*is the receiver location and

_{r}*t*the receiver time. This explains the somewhat complicated notation in (8). The physical interpretation of this result is that the perturbations

_{r}*δa*and

*δb*create scattering sources

*Q*and

_{δa}*Q*, which give rise to a distribution

_{δb}*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

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, 10].

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

$$=0\phantom{\rule{.9em}{0ex}}\mathit{in}\phantom{\rule{.2em}{0ex}}\Omega \times {S}^{n-1}\times \left[0,T\right]$$

*with the* ‘*final value condition*’

*and the inhomogeneous outgoing boundary condition*

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

*with*

*Solving* (15)–(17) *is called* ‘*backtransport*’.

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 *x _{r}* and detection time

*t*. 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

_{r}*∂*and the space derivative

_{t}*θ*·∇ 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

*x*during the experiment, i.e.

_{sc}*u*(

*x*,

_{sc}*θ*,

*t*)=0 for all

*θ*∊

*S*

^{n}^{-1}and all

*t*∊[0,

*T*], then the update (

*R′**

*ζ)(*

_{a,b}*x*) in (19), (20) will be zero at this location

*x*. 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 (

_{sc}*δa*(

*x*),

_{sc}*δb*(

*x*)) 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

_{sc}*x*in the medium. We will discuss these sensitivity functions for the linear transport equation in the following section.

_{sc}## 3 Transport sensitivity functions

The following result provides a formal link between *backtransport*, as it is derived above, and the general idea of *backprojection*, which is well-known in x-ray tomography [17]. A short derivation of this result is given in the appendix.

Transport sensitivity functions. *For a given source (5), we can find functions* ψ*a*(*x _{r}*,

*t*;

_{r}*x*) and ψ

_{sc}*(*

_{b}*x*,

_{r}*t*;

_{r}*x*)

_{sc}*with the following properties*.

1.)‘*Projection*’ *Step*: *The action of the linearized residual operator*
* R′ _{a,b}*

*on the perturbation*((

*δa*,

*δb*)

*in parameter space can be described as follows*

$${\int}_{\Omega}{\mathit{dx}}_{\mathit{sc}}{\Psi}_{a}({x}_{r},{t}_{r};{x}_{\mathit{sc}})\mathit{\delta a}\left({x}_{\mathit{sc}}\right)+{\int}_{\Omega}{\mathit{dx}}_{\mathit{sc}}{\Psi}_{b}({x}_{r},{t}_{r};{x}_{\mathit{sc}})\mathit{\delta b}\left({x}_{\mathit{sc}}\right)$$

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

$$\left({\int}_{\partial \Omega}{\mathit{dx}}_{r}{\int}_{\left[0,T\right]}{\mathit{dt}}_{r}{\Psi}_{a}({x}_{r},{t}_{r};{x}_{\mathit{sc}})\varsigma ({x}_{r},{t}_{r}){,\int}_{\partial \Omega}{\mathit{dx}}_{r}{\int}_{\left[0,T\right]}{\mathit{dt}}_{r}{\Psi}_{b}({x}_{r},{t}_{r};{x}_{\mathit{sc}})\varsigma ({x}_{r},{t}_{r})\right)$$

*The functions* ψ* _{a}*(

*x*,

_{r}*t*;

_{r}*x*) and ψ

_{sc}*(*

_{b}*x*,

_{r}*t*;

_{r}*x*)

_{sc}*are called*‘

*transport sensitivity functions*’.

Formula (21) says that a parameter perturbation (*δa*(*x _{sc}*),

*δb*(

*x*)) at a position

_{sc}*x*with a large sensitivity value ψ

_{sc}_{a,b}(

*x*,

_{r}*t*;

_{r}*x*) will have a relatively strong influence on the data measured at the location

_{sc}*x*and at the time

_{r}*t*. Formula (22), on the other hand, says that a residual at a receiver with position

_{r}*x*, which is detected at time

_{r}*t*, will produce relatively large updates of the parameters (

_{r}*a*,

*b*) in the medium at positions

*x*where the corresponding sensitivity functions ψ

_{sc}*(*

_{a,b}*x*,

_{r}*t*;

_{r}*x*) have large values.

_{sc}We mention that in x-ray tomography the action of the adjoint operator of the linear Radon transform to the data is typically called ‘backprojection’. The physical picture in these situations is very similar. In that application, the data are ‘smeared back’ uniformly over those lines in the imaging domain which have contributed to the data. These lines can be considered as parametrized by the source and receiver positions. One of the most popular reconstruction schemes in x-ray tomography (‘filtered backprojection’) is based on this procedure.

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}*x*,

_{r}*t*;

_{r}*x*). 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.

_{sc}## 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

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

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.

An efficient strategy for combining these updates corresponding to single source positions has been derived in [8, 9, 10, 18, 19]. In that approach, a correction (*δa*, *δb*) is applied immediately to the parameter distribution (*a, b*) after being calculated from the data of one source position, which yields an updated best guess (*a*+*δ _{a, b}*+

*δ*). Then, the data corresponding to the next source position are used to find a new correction (

_{b}*δa*,

*δb*) to these parameters by running one forward and one adjoint problem on the corrected guess. This is done iteratively marching from source to source in some cyclic order, until the method converges. Since with each update the reference parameter distribution (

*a, b*) changes, the described inversion method is

*nonlinear*. For a more detailed description of this single-step adjoint field inversion scheme in the application of optical tomography we refer to [8, 9, 10].

We want to mention that the investigation of sensitivity relations in optical tomography is by far not new. It has been done in different forms in [1, 2, 3, 11, 24]. For alternative approaches to solving related inverse problems in optical tomography we refer to [3, 7, 12, 14] and the references therein. A general overview can be found in [3].

## 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, 10]. Our computational domain Ω is a square of the size 6×6 cm^{2} 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, 10].

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, 15, 20, 21]. 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.

In Figure 2 we have displayed some typical data *M _{a,b}*(

*x*,

_{r}*t*), defined in (6), for these four experiments. In the first image, a delta-like laser source (5) is positioned at the boundary point

_{r}*P*1(which is marked in Figure 1 in the upper left image) with coordinates (40, 0). The receiver is situated at the boundary point

*P*3 (see again Figure 1) and has the coordinates (60, 120). In the second image, the source is positioned at the boundary point

*P*2 with coordinates (120,40), and the receiver is again situated at the boundary point

*P*3. 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, 21] which should provide an interesting alternative to the pure transport model used here for the modelling of scattering media with voids.

In the Figures 3–9 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.

Figure 3 displays -ψ* _{a}*(

*x*,

_{r}*t*;

_{r}*x*) where the source is situated at position

_{sc}*P*1 and the receiver location

*x*is at

_{r}*P*3. The receiving time is

*t*=10 s, which corresponds to the time step

_{r}*T*1=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 *t _{r}*=24 s, which corresponds to time step

*T*2=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

*P*2 instead of

*P*1. 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*(*x _{r}*,

*t*;

_{r}*x*) 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

_{sc}*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.

In Figure 8 we have displayed ψ* _{b}*(

*x*,

_{r}*t*;

_{r}*x*) 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.

_{sc}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, 23], 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, 22].

## 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*[*x _{sc}*,

*θ*,

_{sc}*t*](

_{sc}*x*,

*θ*,

*t*) as solutions of the problem

$$=\delta (x-{x}_{\mathit{sc}})\delta (\theta -{\theta}_{\mathit{sc}})\delta (t-{t}_{\mathit{sc}})$$

in Ω×*S ^{n}*

^{-1}×[0,

*T*] with the homogeneous initial condition

and the absorbing boundary condition

The index sc indicates an arbitrary point in the scattering region. The *adjoint transport Green functions*
*G*̂[*x _{r}*,

*θ*,

_{r}*t*](

_{r}*x*,

*θ*,

*t*) are defined accordingly by

$$=\delta (x-{x}_{r})\delta (\theta -{\theta}_{r})\delta (t-{t}_{r})$$

in Ω×*S ^{n}*

^{-1}×[0,

*T*] with the ‘final value condition’

and the outgoing boundary condition

The functions Φ* _{a}*[

*x*,

_{r}*θ*,

_{r}*t*] are defined by

_{r}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

$$\times G[{x}_{s},{\theta}_{s},0]({x}_{\mathit{sc}},{\theta}_{\mathit{sc}},{t}_{\mathit{sc}})G[{x}_{\mathit{sc}},{\theta}_{\mathit{sc}},{t}_{\mathit{sc}}]({x}_{r},{\theta}_{r},{t}_{r})\mathit{\delta a}\left({x}_{\mathit{sc}}\right)$$

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

and defining the absorption sensitivity functions ψ* _{a}*(

*x*,

_{r}*t*;

_{r}*x*) by

_{sc}we can write (32) in the shorter form

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}*[

*x*,

_{r}*θ*,

_{r}*t*] by

_{r}$$\left[G[{x}_{s},{\theta}_{s},0]({x}_{\mathit{sc}},{\theta}_{\mathit{sc}},{t}_{\mathit{sc}})-{\int}_{{S}^{n-1}}\eta ({\theta}_{\mathit{sc}}\xb7{\theta}^{\prime})G[{x}_{s},{\theta}_{s},0]({x}_{\mathit{sc}},{\theta}^{\prime},{t}_{\mathit{sc}}){\mathit{d\theta}}^{\prime}\right].$$

Introducing the scattering sensitivity functions ψ* _{b}*(

*x*,

_{r}*t*;

_{r}*x*) as

_{sc}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].

Practically, we only have to solve one forward transport problem (1)–(3) and one adjoint transport problem (15)–(17) with ζ,(*x*,*t*)=δ(*x*-*x _{r}*)

*δ*(

*t*-

*t*) in order to calculate the sensitivity functions ψ

_{r}*(*

_{a}*x*,

_{r}*t*;

_{r}*x*) and ψ

_{sc}*(*

_{b}*x*,

_{r}*t*;

_{r}*x*) corresponding to a given receiver position

_{sc}*x*and a receiver time

_{r}*t*. The sensitivity functions are then automatically given by ψ

_{r}*(*

_{a}*x*,

_{r}*t*;

_{r}*x*)=Δ

_{sc}*(*

_{a}*x*) and ψ

_{sc}*(*

_{b}*x*,

_{r}*t*;

_{r}*x*)=Δ

_{sc}*(*

_{b}*x*) with Δ

_{sc}*and Δ*

_{a}*defined in (19) and (20).*

_{b}## 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]