## Abstract

We present an investigation of the effect of a 3D non-scattering gap region on image reconstruction in diffuse optical tomography. The void gap is modelled by the Radiosity-Diffusion method and the inverse problem is solved using the adjoint field method. The case of a sphere with concentric spherical gap is used as an example.

©2000 Optical Society of America

## 1 Introduction

Optical Tomography is an example of a non-linear illposed inverse problem. A complete treatment of image reconstruction therefore requires either implicitly or explicitly the optimization of an objective function with respect to the parameters of a model[1]. The most commonly used model is the diffusion approximation which is the lowest order non-trivial approximation to the more correct Radiative Transfer Equation (RTE)[1, 2]. The reason for its success is that tissues of interest in clinical applications, such as the human breast and peripheral muscles, are of a scale such that detected photons have undergone a high enough number of scattering events that the dependence of fluence on direction is only linear.

However, one major application interest for optical tomography is its use for imaging of the brain. Within the head the regions of cerebrospinal fluid (CSF) around the brain and in the ventricles are non-scattering while still absorbing. The presence of CSF prevents the accurate modeling of photon propagation within the regions of interest when using the diffusion approximation and leads naturally to the attempt to develop a modelling scheme based on the RTE[3, 4, 5, 6]. The principle difficulty then is that accurate handling of the non-scattering regions implies the use of a very computationally expensive numerical method.

One proposed alternative to the full RTE is the *Radiosity-Diffusion* model which assumes diffusive regions coupled by non-scattering voids [7, 8]. Propagation in the void is handled via geometrical optics and the interface requires the development of *non-local* boundary conditions [9]. The advantage is the greatly increased computational efficiency whilst maintaining good accuracy. This efficient model has allowed the investigation of the effect of voids on the accuracy of image reconstruction in 2D [10]. However optical tomography is in reality a 3D problem [11]. In this paper we report the development of the radiosity-diffusion model in 3D and present the first assessment of optical tomography in a 3D geometry involving voids.

## 2 The Radiosity-Diffusion Model

Let Ω represent a domain consisting of *R* diffusing regions {Ω_{1}, Ω_{2}, … Ω
_{R}
} and *V* void regions {Ξ_{1}, Ξ_{2},… Ξ
_{V}
}. Let Ω
^{d}
=${\cup}_{i=1}^{R}$Ω
_{i}
be the union of all diffusing regions and Ξ
^{d}
=${\cup}_{i=1}^{V}$ Ξ
_{i}
be the union of all void regions. Thus Ω=Ω
^{d}
∪ Ξ
^{d}
. Each region has an outer boundary ${\mathit{\partial}}_{i}^{+}$
and a non-negative number of inner boundaries ${\mathit{\partial}}_{i\mathit{,}k}^{-}$. We define ∂Ω
_{i}
=∂${\mathrm{\Omega}}_{i}^{+}$
∪(∪
_{k}
∂${\mathrm{\Omega}}_{i\mathit{,}k}^{-}$
) for diffusing regions, with similarly ∂Ξ
_{i}
=∂${\mathrm{\Xi}}_{i}^{+}$
∪(∪
_{k}
∂${\mathrm{\Xi}}_{i\mathit{,}k}^{-}$
) for voids. The outer boundary ∂${\mathrm{\Omega}}_{1}^{+}$ will also be denoted by *∂*Ω.

We consider a frequency-domain system so that the photon density Φ(*r*; *ω*) at modulation frequency *ω* satisfies the homogeneous equation

where *κ*=1/(3(*µ*
_{a}+*µ′*
_{s})) is the diffusion coefficient defined in terms of *µ*
_{a}(*r*) and *µ′*
_{s}(*r*), the spatially varying absorption and reduced scattering coefficients respectively and with local inhomogeneous boundary conditions

where *η* is the incoming flux modelled as a Neumann source [13], and *v* is the surface normal pointing into the void region, and non-local boundary conditions

$$\frac{\mathrm{exp}\left[-\left({\mu}_{a}+i\frac{\omega}{c}\right)\mid \mathit{m}-{\mathit{m}}^{\prime}\mid \right]}{{\mid \mathit{m}-{\mathit{m}}^{\prime}\mid}^{2}}d{\mathit{m}}^{\prime}\phantom{\rule{.9em}{0ex}}\mathit{m},{\mathit{m}}^{\prime}\u220a\partial \Xi $$

$$\mathrm{cos}\theta =\hat{\mathit{\nu}}\left(\mathit{m}\right)\xb7\frac{\mathit{m}-{\mathit{m}}^{\prime}}{\mid \mathit{m}-{\mathit{m}}^{\prime}\mid}\phantom{\rule{.9em}{0ex}},\mathrm{cos}{\theta}^{\prime}=\hat{\mathit{\nu}}\left({\mathit{m}}^{\prime}\right)\xb7\frac{\mathit{m}-{\mathit{m}}^{\prime}}{\mid \mathit{m}-{\mathit{m}}^{\prime}\mid}$$

and *h*(*m*, *m′*) is unity if *m*, *m′* are in line of sight across the void, and zero otherwise. A more exact boundary condition replaces the term Φ(*m′*;*ω*)/2*A* in eq(3)by

$\left(\Phi ({m}^{\prime};\omega )-2\kappa \left({m}^{\prime}\right)\frac{1+{R}^{\left(1\right)}}{1-{R}^{\left(0\right)}}\frac{\partial \Phi ({m}^{\prime};\omega )}{\partial \nu}\right)(1-{\mid R\left(\theta \right)\mid}^{2})$

where *R*
^{(0)}, *R*
^{(1)} are derived from the Fresnel coefficients taken over local coordinates

*R*
^{(0)}=2${\int}_{\mathrm{\pi}/2}^{\mathrm{\pi}}$
*R*(*ϑ*)sin*ϑ*cosϑdϑ, *R*
^{(1)}=3${\int}_{\mathrm{\pi}/2}^{\mathrm{\pi}}$
*R*(ϑ)sin ϑcos^{2}ϑdϑ

See [12] for a detailed derivation. These two forms were compared in [9] and found to give comparable results. The constant for the domain boundary in eq(2) and (3) is given by *A*=(1-*R*
^{(1)})/(1-*R*
^{(0)}). For each Neumann source term *η* the measureable is

## 3 Inverse Problem

In the adjoint field approach [14] we assume a finite number of input fluxes {*η*_{j}
; *j*=1,…, *S*}, and given data {*g*_{j}
;*j*=1,…, *S*}, and look for the minimisation of the norm

A descent direction for minimising *C* with respect to the data from the *j*^{th}
source is given by

with ψ
_{j}
the solution to the adjoint diffusion equation

with boundary conditions

$${\Psi}_{\mathit{j}}(m;\mathit{\omega})+2A\kappa \frac{\partial {\Psi}_{j}(m;\omega )}{\partial \mathit{\nu}}=\frac{1}{\pi}{\int}_{\partial \Xi}\mathrm{cos}\theta \mathrm{cos}\theta \text{'}\frac{{\Psi}_{j}({m}^{\prime};\omega )}{2A}h(m,{m}^{\prime})x$$

$$\frac{\mathrm{exp}[-({\mu}_{a}-i\frac{\omega}{c})\mid m-{m}^{\prime}\mid ]}{{\mid \mathit{m}-{\mathit{m}}^{\prime}\mid}^{2}}d{m}^{\prime}\phantom{\rule{.9em}{0ex}}m,{m}^{\prime}\u220a\partial \Xi $$

For the optimization scheme we employ the same conjugate-gradient scheme as used in [10]. However, adequate treatment of the 3D problem requires some special attention, which we describe in the next section.

## 4 Implementation in 3D

Our numerical method utilises the Finite Element Method (FEM) to allow application to general domain boundaries and parameter distributions. In the results presented here we use linear tetrahedral elements generated using the public domain package Netgen [15], which implies that the fluence Φ is expressed as a linear combination of piecewise linear shape functions {*u*_{k}
; *k*=1… *D*} with value one at nodal points *N*_{k}
and zero at all other nodes. For simplicity we consider the case of concentric spherical surfaces which are defined parametrically. The discretisation of eq(1) leads to sparse symmetric matrices, whilst the implementation of the non-local boundary conditions eq(3) leads to a dense *coupling matrix* whose order, *D*_{B}
, is the number of nodes in the void boundary. In 3D this number is proportional to the surface area of the void surfaces which is the principle extra overhead to the computational requirements of the foward solver. In the previous 2D treatment [8] we calculated the visibility function by a full search of all void nodes which becomes computationally expensive in 3D. For the case considered here we can determine the visibility function *h*(*m*, *m′*) as well as the surface normal functions *$\widehat{\nu}$* directly from the parametric representation of the surfaces which greatly simplifies the setup time of the forward solver. The function

extracted from eq(3), is called the *form factor* corresponding to its use in Computer Graphics [16]. In our current implementation we approximate this function by a bilinear expansion over the FEM shape functions on the surface of the void [17]. Between two element faces *τ*_{α}
, *τ*_{α′}
; ∊ *δ*Ξ containing *N*, *N′* nodes respectively the form factor can be approximated by

where *n*(*k*) maps the local node to the global nodes in the complete mesh.

## 5 Results

In fig.1 we show the geometry being considered. We take a sphere radius 25*mm* with concentric sphere radius *r*_{inner}
and a concentric void gap with outer radius *r*_{outer}
. 32 sources and detectors are arranged in three rings. We used background parameters *µ*
_{a}=0.01 *mm*
^{-1}, *µ′*
_{s}=1 *mm*
^{-1} with *µ*
_{a}=0.005 in the void region. For the sources and detectors we use a cosine weighted patch in the parameters of the surface representation. For the reconstruction basis we use (20×20×20) tricubic interpolated voxels.

In fig.2 we show the surface data generated by a single source for the case of a solid sphere and for three gaps defined by *r*_{inner}
=17,16,15*mm*, and in each case *r*_{outer}
=20*mm*. The difference between each void case and the solid sphere emphasises the increased light intensity and decreased mean time in the gap region.

In the animation attached to fig.3 we show the sensitivity functions in absorption for intensity (*y*(*ω*=0)) and mean time $\frac{\partial \mathrm{ln}\mathit{y}}{\partial \omega}{\mid}_{\omega =0}$ measurement types, in the case of a 3 mm wide gap, for 16 different detector locations around the equator of the sphere. The sensitivity functions are derived by using eq(6) with ψ
_{j}
an adjoint Green’s function obtained using a *δ*-function for the local boundary condition in eq(8); see [1, 18] for a more detailed derivation, and the related paper [6]. For intensity measurements the sensitivity is mainly confined to the regions around the source and detector sites in the outer ring, whereas for mean time the gap has less of an adverse effect, providing higher sensitivity in the inner region. We use mean time data only in the following for consistency with [10].

In fig.4 we show the reconstruction of a blob in the 3*mm* gap. The blob had a radius of 3*mm* and an absorption of 0.02*mm*
^{-1} with its center placed at position (12,0,0). The images shown are at iteration 40 of the conjugate gradient scheme. The reconstruction shows good localisation despite quite a large degree of broadening.

## 6 Discussion and conclusions

The correct treatment of a non-scattering region in diffusing media is an important requirement for application of Optical Tomography to the brain. In this paper we showed initial results of a 3D extension to the Radiosity-Diffusion model that provides one approach to this problem. Even though the cases considered here are simple, they indicate that, if the location of the void is known, the signal detected on the domain surface still carries enough information to allow a reconstruction. A more complete treatment should consider data types in addition to just the mean time as used here, as well as simultaneous reconstruction of absorption and scattering. Furthermore, one of the interesting results found in [10] was the improvement in reconstruction quality when the boundary of the void gap is not smooth. Investigation of this case requires a more sophisticated determination of form factors which is under development.

In the case where the void location is not known, we are investigating the use of a boundary recovery scheme related to the method reported in [19]. The application of this scheme to the Radiosity-Diffusion model requires the determination of differentiability of the forward operator at the diffuse/non-diffuse interface.

## Acknowledgments

We thank the Wellcome Trust and the EPRSC. J.Riley thanks UCL Graduate School for his PhD scholarship. J.Ripoll acknowledges a grant from Ministerio de Educación y Cultura. N.Nieto-Vesperinas thanks DGICYT and Fundación Ramón Areces. We thank F. Natterer and O. Dorn for discussions.

## References and links

**1. **S. R. Arridge, “Optical Tomography in Medical Imaging,” Inverse Problems **15**, R41–R93 (1999). [CrossRef]

**2. **S. R. Arridge and J. C. Hebden, “Optical Imaging in Medicine: II. Modelling and Reconstruction,” Phys. Med. Biol. **42**, 841–853 (1997). [CrossRef] [PubMed]

**3. **A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of Finite-Difference Transport and Diffusion Calculations for Photon Migration in Homogeneous and Hetergeneous Tissue,” Phys. Med. Biol. **43**, 1285–1302 (1998). [CrossRef] [PubMed]

**4. **O. Dorn, “A Transport-BackTransport Method for Optical Tomography,” Inverse Problems **14**, 1107–1130 (1998). [CrossRef]

**5. **A. D. Klose and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. **26**1698–1707 (1999). [CrossRef] [PubMed]

**6. **O. Dorn, “Scattering and absorption transport sensitivity functions for optical tomography,” Opt. Express This issue (2000).

**7. **M. Firbank, S. R. Arridge, M. Schweiger, and D. T. Delpy, “An investigation of light transport through scattering bodies with non-scattering regions,” Phys. Med. Biol. **41**, 767–783 (1996). [CrossRef] [PubMed]

**8. **S. R. Arridge, H. Dehghani, 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,” Med. Phys. **27**, 252–264 (2000). [CrossRef] [PubMed]

**9. **J. Ripoll, S. R. Arridge, H. Dehghani, and M. Nieto-Vesperinas, “Boundary conditions for light propagation in diffusive media with nonscattering regions,” J. Opt. Soc. Am. A **17**, 1671–1681 (2000). [CrossRef]

**10. **H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical Tomography in the Presence of Void Regions,” J. Opt. Soc. Am. A **17**, 1659–1670 (2000). [CrossRef]

**11. **M. Schweiger and S. R. Arridge, “Comparison of 2D and 3D reconstruction algorithms in Optical Tomography,” Appl. Opt. **37**, 7419–7428 (1998). [CrossRef]

**12. **J. Ripoll, Ph.D. thesis, University Autónoma of Madrid, 2000.

**13. **S. R. Arridge and M. Schweiger, “The Finite Element Model for the Propagation of Light in Scattering Media: Boundary and Source Conditions,” Med. Phys. **22**, 1779–1792, (1995). [CrossRef] [PubMed]

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

**15. **J. Schoberl, “NetGen”, http://www.sfb013.uni-linz.ac.at/joachim/netgen/

**16. **M. F. Cohen and J. R. Wallace, *Radiosity and Realistic Image Synthesis* (Academic, London, 1993).

**17. **H. R. Zatz, Master’s thesis, Cornell University, 1993.

**18. **S. R. Arridge and M. Schweiger, “Photon Measurement Density Functions. Part 2: Finite Element Calculations,” Appl. Opt. **34**, 8026–8037 (1995). [CrossRef] [PubMed]

**19. **V. Kolehmainen, M. Vaukhonen, J. P. Kaipio, and S. R. Arridge, “Recovery of piecewise constant coefficients in optical diffusion tomography,” Opt. Express This issue (2000).