1 The transport equation in optical tomography

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

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=2,3) with smooth boundary ∂Ω. S ^n-1 denotes the set of direction vectors θ which have unit length in ℝⁿ. 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_s by 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_a≈1.0–10.0 cm, 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 ·{\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 θ_s. The symbols δ(x-x_s), δ(t) and δ(θ-θ_s) denote Dirac delta functions in the corresponding variables. If the source is positioned at the boundary, x_s∊∂Ω, we assume that v(x_s)·θ_s<0. In our numerical experiments in section 5 we will use v(x_s)·θ_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].

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_r denote the receiver location and receiving time, respectively. The symbol 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

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_sc)), where the index sc stands for ‘scatter’, the value of R′_a,b(δa, δb) is a function in the variables x_r and t_r, where x_r is the receiver location and t_r 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

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

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_r. 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 x_sc during the experiment, i.e. 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_sc. 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(x_sc), δb(x_sc)) 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 x_sc in the medium. We will discuss these sensitivity functions for the linear transport equation in the following section.

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_sc) and ψ_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

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

The functions ψ_a(x_r,t_r;x_sc) and ψ_b(x_r,t_r;x_sc) are called ‘transport sensitivity functions’.

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

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_sc). 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

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}+δ_b). Then, the data corresponding to the next source position are used to find a new correction (δ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].

 

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.

Download Full Size | PDF

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² 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 ¹, 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.

 

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.

Download Full Size | PDF

In Figure 2 we have displayed some typical data M_a,b(x_r, t_r), 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, 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.

 

Fig. 3. Transport absorption sensitivity functions -ψ_a(x_r, t_r; x_sc) for a source located at P1 and a receiver at P3. The receiving time is t_r=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).

Download Full Size | PDF

Figure 3 displays -ψ_a(x_r, t_r; x_sc) where the source is situated at position P1 and the receiver location x_r is at P3. The receiving time is t_r=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 t_r=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(x_r, t_r; x_sc) 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(x_r, t_r; x_sc) for a source located at P1 and a receiver at P3. The receiving time is t_r=24 s, which corresponds to the time step T2=120. Top left: experiment 1 (N=3.0×10⁴); Top right: experiment 2 (N=5.4×10⁴); Bottom left: experiment 3 (N=7.3×10⁴); Bottom right: experiment 4 (N=4.1×10⁴).

Download Full Size | PDF

 

Fig. 5. Transport absorption sensitivity functions -ψ_a(x_r, t_r; x_sc) for a source located at P2 and a receiver at P3. The receiving time is t_r=10 s, which corresponds to the time step T1=50. Top left: experiment 1 (N=5.6×10³); Top right: experiment 2 (N=8.5×10⁴); Bottom left: experiment 3 (N=1.0×10⁵); Bottom right: experiment 4 (N=1.3×10⁴).

Download Full Size | PDF

 

Fig. 6. Transport absorption sensitivity functions -ψ_a(x_r, t_r; x_sc) for a source located at P2 and a receiver at P3. The receiving time is t_r=24 s, which corresponds to the time step T2=120. Top left: experiment 1(N=9.8×10⁴); Top right: experiment 2(N=1.7×10⁵); Bottom left: experiment 3(N=2.3×10⁵); Bottom right: experiment 4 (N=1.4×10⁵).

Download Full Size | PDF

In Figure 8 we have displayed ψ_b(x_r, t_r; x_sc) 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, 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

(25)$$\frac{\partial G}{\partial t}+\theta ·\nabla G(x,\theta ,t)+\left(a\left(x\right)+b\left(x\right)\right)G(x,\theta ,t)-b\left(x\right){\int }_{S^{n-1}}\eta (\theta ·{\theta }^{\prime })G(x,{\theta }^{\prime },t){\mathit{d\theta }}^{\prime }$$
$$=\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

(26)$$G(x,\theta ,0)=0\mathrm{in}\Omega \times s^{n-1},$$

and the absorbing boundary condition

(27)$$G(x,\theta ,t)=0\mathrm{on}{\Gamma }_{-}.$$

 

Fig. 7. Transport scattering sensitivity functions ψ_b(x_r, t_r; x_sc) for a source located at P1 and a receiver at P3. The receiving time is t_r=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).

Download Full Size | PDF

 

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(x_r, t_r; x_sc) for a source located at P1 and a receiver at P3. The receiving time step is T1=50.

Download Full Size | PDF

 

Fig. 9. Transport scattering sensitivity functions ψb(x_r, t_r; x_sc) for a source located at P1 and a receiver at P3. The receiving time is t_r=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).

Download Full Size | PDF

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

(28)$$\frac{\partial \hat{G}}{\partial t}-\theta ·\nabla \hat{G}(x,\theta ,t)+\left(a\left(x\right)+b\left(x\right)\right)\hat{G}(x,\theta ,t)-b\left(x\right){\int }_{S^{n-1}}\eta (\theta ·{\theta }^{\prime })\hat{G}(x,{\theta }^{\prime },t){\mathit{d\theta }}^{\prime }$$
$$=\delta (x-x_r)\delta (\theta -{\theta }_r)\delta (t-t_r)$$

in Ω×S^{n -1}×[0,T] with the ‘final value condition’

(29)$$\hat{G}(x,\theta ,T)=0\mathrm{in}\Omega \times S^{n-1},$$

and the outgoing boundary condition

(30)$$\hat{G}(x,\theta ,t)=0\mathrm{on}{\Gamma }_{+}.$$

The functions Φ_a[x_r, θ_r, t_r] are defined by

(31)$${\Phi }_a[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}})=\hat{G}[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathrm{sc}})G[x_s,{\theta }_s,0](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}}).$$

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

(32)$$\left(R_{a,b}^{\prime }(\mathit{\delta a},0)\right)(x_r,t_r)=-{\int }_{S_{+}^{n-1}}{\mathit{d\theta }}_r{\int }_{\Omega }{\mathit{dx}}_{\mathit{sc}}{\int }_{S^{n-1}}{\mathit{d\theta }}_{\mathit{sc}}{\int }_{\left[0,T\right]}{\mathit{dt}}_{\mathit{sc}}\nu \left(x_r\right)·{\theta }_r$$
$$\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

(33)$$G[x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}}](x_r,{\theta }_r,t_r)=\hat{G}[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},{\mathit{t}}_{\mathit{sc}}),$$

and defining the absorption sensitivity functions ψ_a(x_r,t_r;x_sc) by

(34)$${\Psi }_a(x_r,t_r;x_{\mathit{sc}})=-{\int }_{S_{+}^{\mathit{n}-1}}{\mathit{d\theta }}_r\nu \left(x_r\right)·{\theta }_r{\int }_{\left[0,T\right]}{\mathit{dt}}_{\mathit{sc}}{\int }_{S^{\mathit{n}-1}}{\mathit{d\theta }}_{\mathit{sc}}{\Phi }_a[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}}),$$

we can write (32) in the shorter form

(35)$$\left(R_{a,b}^{\prime }(\mathit{\delta a},0)\right)(x_r,t_r)={\int }_{\Omega }{\mathit{dx}}_{\mathit{sc}}{\Psi }_a(x_r,t_r;x_{\mathit{sc}})\mathit{\delta a}\left(x_{\mathit{sc}}\right),$$

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_r] by

(36)$${\Phi }_b[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}})=\hat{G}[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}})\times$$
$$\left[G[x_s,{\theta }_s,0](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{sc}})-{\int }_{S^{n-1}}\eta ({\theta }_{\mathit{sc}}·{\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_sc) as

(37)$${\Psi }_b(x_r,t_r,x_{\mathit{sc}})={\int }_{S_{+}^{n-1}}{\mathit{d\theta }}_r\mathit{\nu }\left(x_r\right).{\mathit{\theta }}_r{\int }_{\left[0,T\right]}{\mathit{dt}}_{\mathit{sc}}{\int }_{S^{n-1}}{\mathit{d\theta }}_{\mathit{sc}}{\Phi }_b[x_r,{\theta }_r,t_r](x_{\mathit{sc}},{\theta }_{\mathit{sc}},t_{\mathit{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_r) in order to calculate the sensitivity functions ψ_a(x_r, t_r; x_sc) and ψ_b(x_r, t_r; x_sc) corresponding to a given receiver position x_r and a receiver time t_r. The sensitivity functions are then automatically given by ψ_a(x_r,t_r; x_sc)=Δ_a(x_sc) and ψ_b(x_r, t_r; x_sc)=Δ_b(x_sc) 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]