## Abstract

Diffuse optical tomography (DOT) aims at recovering three-dimensional images of absorption and scattering parameters inside diffusive body based on small number of transmission measurements at the boundary of the body. This image reconstruction problem is known to be an ill-posed inverse problem, which requires use of prior information for successful reconstruction. We present a shape based method for DOT, where we assume *a priori* that the unknown body consist of disjoint subdomains with different optical properties. We utilize spherical harmonics expansion to parameterize the reconstruction problem with respect to the subdomain boundaries, and introduce a finite element (FEM) based algorithm that uses a novel 3D mesh subdivision technique to describe the mapping from spherical harmonics coefficients to the 3D absorption and scattering distributions inside a unstructured volumetric FEM mesh. We evaluate the shape based method by reconstructing experimental DOT data, from a cylindrical phantom with one inclusion with high absorption and one with high scattering. The reconstruction was monitored, and we found a 87% reduction in the Hausdorff measure between targets and reconstructed inclusions, 96% success in recovering the location of the centers of the inclusions and 87% success in average in the recovery for the volumes.

©2009 Optical Society of America

## 1. Introduction

In this paper we present a shape based reconstruction method for Diffuse Optical tomography (DOT). DOT [1], seeks the recovery of absorption and scattering parameters of biological tissues, given measurements of transmitted light through tissue of several centimeters in thickness. There are several physiologically interesting parameters which can be derived from the knowledge of the absorption and scattering coefficients of light in tissue. These include tissue oxygenation, blood volume and blood oxygenation [2, 3]. Primary applications are the detection and classification of tumorous tissue in the breast, monitoring of the oxygenation level in infant brain tissue, and functional brain activation studies.

Image reconstruction in DOT is a non linear inverse problem, which is highly unstable with respect to measurement noise. Traditionally, this image reconstruction problem is solved using a high dimensional voxel based parameterization for the absorption and scattering parameters. The success of the reconstruction depends on the inclusion of additional (*a priori*) information to stabilize the reconstruction. Most commonly the prior information is included into the problem either in the form of regularisation techniques, which imply the specification of a prior probability distribution of voxel values based on geometric, local smoothness, or statistical asumptions[4, 5, 7, 8, 15], by specification of tissue classifications and projection into a subspace formed by regional groupings[6, 10, 11, 12], or a combination of methods[13, 14]. Reviews on the approaches for inclusion of a priori information can be found, for example, in [9, 16, 17].

In this paper, we consider a different approach where the unknown body can be assumed to consist of a small number of subdomains of different tissues with different absorption and scattering values. Based on the assumption that the unknowns in such cases could be defined by the subdomain boundaries, a field of shape based reconstruction techniques has emerged in recent years. There have been many results in literature of either level-set techniques [18, 19, 20, 21, 22, 23, 24], where the boundaries of the subdomains are implicitly modeled by the zero level of level-set functions, or parametric methods [25, 26, 27, 28] where an explicit parameterisation of the boundaries is used. In the case of the shape based reconstructions that use an explicit parameterisation of the boundaries, many different ways to describe the shapes have been used such as spherical harmonics, ellipsoids and spheres in 3D or Fourier curves, and Hermite polynomials, in 2D. The use of shape parameterisation makes the reconstruction problem less ill-posed by providing implicit regularisation due to the reduced dimensionality of the problem. Furthermore, in shape based approaches *a priori* information about the topology, the approximate locations and shapes of the unknown subdomains and their optical properties can be incorporated in the inversion directly and the division of the body to different tissue types is obtained without post process segmentation.

In this paper, we present a three-dimensional (3D) shape based approach for DOT. We parameterize the inverse problem with respect to the spherical harmonics coefficients of the subdomain boundaries. To develop the FEM based forward mapping from the spherical harmonics coefficients to the DOT data, we generalize the 2D element division technique presented in [26] into a 3D element subdivision technique for the mapping of the spherical harmonics coefficients to the absorption and scattering distributions inside an unstructured volumetric FEM mesh. In our previous work [25] we have evaluated a boundary element method (BEM) based 3D shape estimation method for reconstructing simulated DOT data. In this work we evaluate the shape based method by reconstructing experimental data and recover the location and shape of absorption and scattering targets inside a cylindrical phantom.

This paper is organized as follows. In section 2 we present the diffusion equation based model used for the transport of light in DOT. Section 3 describes the discretisation steps for modeling of measurements for given shape parameters describing the boundaries of the subdomains. In section 4 the inverse problem is formulated as finding the spherical harmonic coefficients that minimize the least squares residual between the experimentally measured and the predicted data. The derivation of the Jacobian of the forward model and an iterative method for solving the inverse problem is implemented in the following subsections. Experimental phantom results are given in 5. Finally, conclusions and ideas for future work are given in the last section of this paper.

## 2. Forward model

The most common model for OT in a highly scattering medium is the diffusion approximation to the radiative transfer equation [1]. In a domain Ω bounded by boundary *∂*Ω we consider the diffusion equation in the frequency case :

In (1), *D* and *µ*
_{a} are the diffusion and absorption coefficients, *ω* is the source modulation frequency, Φ is the photon density and c the speed of light in the medium. The boundary source distribution *q*(**r**;*ω*) represents the number of photons at a boundary position **r**; The diffusion coefficient is given by

Where *µ*′_{s} denotes the reduced scattering coefficient. The appropriate boundary condition for our problem is of the Robin type

where *ζ* models the refractive index mismatch at the boundary *∂*Ω and *ν* is the outward surface normal at **r**.

In the remainder of the paper we refer to the space of parameters **p**={*µ*
_{a},*D*} as *𝓟*. The space of fields Φ will be denoted as *𝓕* and the space of data y will be denoted as *𝓩*. The data y will be calculated by the use of a measurement operator *𝓕* : *𝓕*→*𝓩* in the form:

The operator typically calculates a vector of weighted integrals of the (complex-valued) fields Φ along a set of surface elements at *∂*Ω where the detectors are located (one integral per detector). For our analysis it is sufficient to assume that this operator is linear and the adjoint measurement operator *M** : *𝓩*→*𝓕* is also well defined. The adjoint operator will be used for defining artificial ‘adjoint sources’ for the calculation of the sensitivity functions in section 4.

Equations (1),(3) and (4) can also be formally described as an operator *A* mapping from the parameter space to the measurement space

## 3. Discretization with FEM

There have been many approaches, such as finite differences, finite elements and boundary elements (BEM), for the solution of the diffusion equation model in arbitrary shaped domains Ω. In this work we have adopted the finite elements method (FEM) to the application of a three dimensional shape based problem. This choice was motivated by the work [26]. In addition in a FEM problem we can easily combine piecewise constant regions together with inhomogeneous background, and the size of the system matrix does not depend on the number of disjoint regions presented in the problem as in the case of BEM, [25].

In FEM the domain Ω is divided into N disjoint elements Ω=∪^{N}
_{n=1}
*H _{n}* joined at

*i*=1…

*h*vertex nodes

*T*. In this paper we present the method in terms of tetrahedral elements, although it could easily be generalized for other types of elements. Given the set of

_{i}*T*nodes and a set of associated basis functions υ

_{i}_{i}(

**r**), the solution for the field Φ defined in Ω, can be approximated by a piecewise continuous function Φ(

**r**)=∑

^{h}

_{i=1}

*ϕ*υ

_{i}_{i}(

**r**), where Φ={

*ϕ*} is the h-dimensional vector of basis coefficients and υi are linear basis functions with limited support over the elements that have the node

_{i}*T*as one of their vertices. From previous work [29, 34], this basis representation transforms the continuous problem of equation (1) and (3) into the linear system of the FEM framework.

_{i}where *α* comes from the boundary conditions (3), and the system matrices K,C,E,B are given by:

$${C}_{ij}={\int}_{\Omega}{\mu}_{a}\left(\mathbf{r}\right){\upsilon}_{i}\left(\mathbf{r}\right){\upsilon}_{j}\left(\mathbf{r}\right)d\mathbf{r}$$

$${E}_{ij}={\int}_{\partial \Omega}{\upsilon}_{i}\left(\mathbf{r}\right){\upsilon}_{j}\left(\mathbf{r}\right)\mathrm{dS}$$

$${B}_{ij}=\genfrac{}{}{0.1ex}{}{1}{c}{\int}_{\Omega}{\upsilon}_{i}\left(\mathbf{r}\right){\upsilon}_{j}\left(\mathbf{r}\right)d\mathbf{r}$$

Equation (6) is formally solved by matrix inversion:

In the shape based-method we assume that the domain Ω contains *L* simply connected subdomains {Ω_{ℓ}}, *ℓ*=1,…,*L*, which are bounded by closed surfaces {Γ_{ℓ}} and have constant optical parameters *µ*
_{a} and *D* (the remaining “background” domain is denoted by Ω_{0}=Ω\∪^{L}
_{ℓ=1}Ω_{ℓ}). Since this is a three dimensional application we choose a real valued spherical harmonics parametric representation for the boundaries of the regions {Ω_{ℓ}} as introduced in [25]. The surface locations r|Γ_{ℓ}, of a surface Γ_{ℓ}, are given by the representation

with expansion coefficients {*C ^{m}_{l}*}. Here, the basis functions

*Y*̃

*(ϑ,φ) are defined as*

^{m}_{l}where *Y ^{m}_{l}* are the (complex-valued) spherical harmonics functions, and

*W*is the maximum degree of spherical harmonics used for the particular representation. For simplicity we introduce the notation

such that the new index *k* ranges over *k*=1,⋯3*W*
^{2}. {*γ ^{ℓ}_{k}*} describes the finite set of spherical harmonics coefficients for the surface Γ

_{ℓ}up to degree

*W*.

Having divided the domain Ω into subregions Ω_{ℓ} with piecewise constant parameters, *µ*
_{a} and *D* can be expressed as

where *χ _{ℓ}* is the characteristic function for a region Ω

_{ℓ}. Substituting in equation (7), and since only the matrices K and C are influenced by the parameters inside the domain, we get:

$${C}_{\mathrm{ij}}=\sum _{\ell =0}^{L}{\int}_{\mathrm{supp}\left({\upsilon}_{i}{\upsilon}_{j}\right)\cap {\Omega}_{\ell}}{\mu}_{a,\ell}\left(\mathbf{r}\right){\upsilon}_{i}\left(\mathbf{r}\right)\xb7{\upsilon}_{j}\left(\mathbf{r}\right)d\mathbf{r}$$

Where supp(*υ _{i}υ_{j}*) represents the part of the domain where both the basis functions are non-zero, that is, the elements that contain both the

*T*and

_{i}*T*nodes.

_{j}The goal now is to express the forward operator as a mapping from the parameter space *𝓟*: {*γ,µ*
_{a},*D*} to the space of data *𝓩*. In the next subsections, we will present the methodology used to convert the FEM algorithm to implement this mapping in five stages.

#### 3.1. Classification of FEM mesh nodes

The first step in our computation is to identify which of the nodes *T _{i}*=(

*x*), of the finite element discretisation are inside a given closed surface Γ

_{i},y_{i},z_{i}_{ℓ}(

*θ,ϕ*). To reduce the computational cost we calculate first the bounding box for the surface Γ

_{ℓ}(

*θ,ϕ*) and work only with the nodes that are inside the bounding box. Using the Jordan surface theorem [33], and the odd-even test we need to identify the intersections of a line starting from a particular node

*T*inside the bounding box with the surface. We discetise the surface Γ

_{i}_{ℓ}(

*θ,ϕ*) using the parametric mapping function, that was initially introduced in [25], to map a triangular mesh originally defined on a sphere to the parametric surface. The surface is now a finely defined triangular approximation Γ

_{ℓ}, with triangles

*G*and nodes

^{ℓ}_{k}*S*. We draw the segment

^{ℓ}_{m}*R*={(

_{i}*x*), (

_{i},y_{i},z_{i}*x*)}, where

_{i},y_{i}, z_{min}*z*is the lower z-coordinate of the volume and check for intersections with each of the triangles

_{min}*G*in the surface approximation Γ

^{ℓ}_{k}_{ℓ}. We calculate the bounding quadrilateral along the

*xy*-plane for each of the triangles in the surface, and check if the projection of the node

*T*along the xy-plane lies inside. For the triangles satisfying this check we further check for intersections with the line segment

_{i}*R*, using the [32] ray-triangle intersection algorithm.

_{i}#### 3.2. Classification of mesh elements

Once the mesh nodes *T _{i}* are classified as inside or outside the region bounded by a parametric surface Γ

_{ℓ}(

*θ,ϕ*), we proceed with the classification of the elements

*H*of the mesh using their respective nodes. We mark an element as

_{n}*outside*if all the nodes are outside the surface, and similarly we mark an element as

*inside*if all nodes are positioned internal to the region. In the case of an element with, one, two or three nodes inside a given region Γ

_{ℓ}we mark it as intercepted. Let

*I*(Γ

_{ℓ}) denote the intersecting elements to the surface Γ

_{ℓ}.

The parameters for both the *outside* and the *inside* elements are constant over the element, so that the elements contributions for the system matrix can be calculated using the equations (7), directly. In the case of an *intercepted* element a more complicated procedure must be performed as we will see in the next subsections.

#### 3.3. Determination of intersections of a surface with the elements of the FEM mesh

In the case of an element *H _{v}*∈

*I*(Γ

_{ℓ}(

*θ,ϕ*)), defined by the nodes {

*T*

^{v}_{1},

*T*

^{v}_{2},

*T*

^{v}_{3},

*T*

^{v}_{4}}, we need to find the intersection points

*k*along the edges of the tetrahedron. Using an element

^{v}*G*∈Γ

_{s}_{ℓ}, with its center the closest to the center of the element

*H*, and using the information from the mesh mapping function we get the spherical coordinates (

_{v}*ϑ*

_{0},

*φ*

_{0}) for the node

*S*

_{0}of the element

*G*. The intersection point

_{s}*k*of the parametrically defined surface with the edge of the tetrahedron between the nodes

_{v}*T*

^{v}_{1}and node

*T*

^{v}_{2}by definition is on that edge, so that:

where we used the cross product of the segments defined by the respective end points. At the same time the intersection point *k _{v}* has to be on the parametrically defined surface. There exists a pair of unknown spherical coordinates ${(\vartheta}_{{k}_{\nu}},{\phi}_{{k}_{\nu}})$ that correspond to this intersection point according to (9)

Solving equation (15) in respect to the spherical coordinates of the unknown intersection point *k _{v}* given by (16), is done using a generic Gauss-Newton non-linear solver with initial values (

*ϑ*

_{0},

*φ*

_{0}). Since those initial values were constructed to be close to final solution the solution is rapid.

#### 3.4. Subdivision of the tetrahedron in subtetrahedra

Once we have found the intersection point along the edges of a *intercepted* tetrahedron we can perform a subdivision of the element into subtetrahedra so that the boundary defined by the spherical harmonics surface Γ_{ℓ}(*θ,ϕ*) is accurately represented in the numerical scheme. For an intercepting element *H _{v}*∈

*I*(Γ

_{ℓ}(

*θ,ϕ*)), defined by the nodes {

*T*

^{v}_{1},

*T*

^{v}_{2},

*T*

^{v}_{3},

*T*

^{v}_{4}} we classify two main cases for the intersection with a surface.

**• Case 1:** One of the nodes, *T ^{v}*

_{(+)}, is on a different side of the parametric surface from the other three {

*T*

^{v}_{(-, i)}}

*i*=1,2,3. In Fig. 1(a) we have denoted this partition of the nodes with the signed distance for each node from the surface : we denote the single node with (+) and the three nodes on the other side of the partition with (-). We should mention that the intersection is treated symmetrically and irrespective of which side of the surface is inside or outside. The node

*T*

^{v}_{(+)}together with the three intersection points {

*k*

^{v}_{1},

*k*

^{v}_{2},

*k*

^{v}_{3}}, calculated according to the previous subsection, form one subtetrahedron. The rest of the element that lies on the negative side of the surface has now the shape of a pentahedron and is split into seven subtetrahedra as shown in Fig. 1(c), where we have exploded the subelements for better comprehension.

**• Case 2** Two of the nodes, *T ^{v}*

_{(+, 1)}

*T*

^{v}_{(+, 2)}are on the positive side of the parametric surface and the other two {

*T*

^{v}_{(-,i)}}

*i*=1,…,2 on the negative. In Fig. 1(b) we denote two nodes with (+) and two with (-). Having found in the previous subsection the four intersection points {

*k*

^{v}_{1},

*k*

^{v}_{2},

*k*

^{v}_{3},

*k*

^{v}_{4}}, of the parametric surface with the edges of the element connecting each of the (+) nodes to the(-) nodes, we divide the tetrahedron into two prismatic solids of five faces each. In sequence the divided substructures are subdivided into three tetrahedra each producing the exploded subtetrahedra of Fig. 1(d).

#### 3.5. Computation of the system matrix

As we have mentioned already the contributions in the system matrix from the elements that belong in wholly to the background or are wholly internal to a particular region are calculated using the integrals (13). For those elements *I*(Γ_{ℓ}(*θ,ϕ*)), that are intersected by the surface of a region we have to calculate the separate integrals ∫∇*υ _{i}*·∇

*υ*

_{j}*d*

**r**and ∫υi·υj dr for each of their subelements. The integrals are mapped from the subtetrahedron to the local element and then evaluated with a Gaussian quadrature in the local element, according to a 3D version of the procedure described in [26].

## 4. Inverse Problem

Given g the vector of experimental measurements, the inverse problem is defined as the minimization of the least squares cost functional,

The shape estimation problem then becomes the one of finding the spherical harmonics parameters that define the boundary configuration that minimizes the functional above, that is

A typical way to minimize such a cost function is a Newton-type method, [30, 35], where we search for a minimum for Ξ(*γ*) by iterations of local linearization and Taylor expansion around the current estimate as

where Λ, is a regularisation matrix.

#### 4.1. Jacobian calculation

The Jacobian matrix, in respect to the shape parameters, is defined as

Looking at (6) and (13) and noting:

we obtain:

which leads to

By applying the discrete version of the measurement operator to both sides of the equation we obtain the Jacobian

where Φ* denotes the adjoint solution. Here, we should note that the exterior boundary *∂*Ω is considered known and the shape parameters do not influence the external boundary. The inclusions are also assumed mutually disjoint and disjoint of the exterior boundary. Therefore the derivative of S in respect to the shape parameters *γ _{k}*are not dependent on the boundary term E or the frequency term B. The elements of the derivatives of S with respect to the shape parameters

*γ*are derived in equations (22)–(32)of [26].

_{k}where (*δµ*
_{a}) and (*δD*) denote the difference of *µ*
_{a} and *D* between inside and outside the boundary Γ. Using (25) on (24) we get the element *J _{sd,k}*, corresponding to the source

*s*and detector

*d*for the shape parameter

*k*, of the Jacobian:

where

are the integrals over surface Γ of the shape coefficient *γ _{k}* that describe that surface.

From the derivation of the Jacobian above we note that integrations over the surfaces of the boundaries Γ_{ℓ} of the regions Ω* _{ℓ}* are necessary. For the calculation of the fields Φ

_{ℓ}and the adjoint fields Φ*

_{ℓ}on the nodes

*S*of the surface Γ

_{ℓ}_{ℓ}we find first the element

*E*of the volume discretisation where the node

_{p}*S*of the surface belongs to, and recall the shape functions of this element.

_{s}where *l* ranges through the number of nodes in one element in the FEM discretisation (4 for a tetrahedra).

#### 4.2. Scalings

To improve the conditioning of the inverse problem we use rescaling of the data [34] : we take the log of the frequency domain complex data and also separate the real and imaginary part, so that:

which reflects to the scaled Jacobian J̃ defined as

## 5. Results

A cylindrical phantom of diameter 69.25 mm and height 110 mm made from epoxy resin was used to acquire data using a frequency-domain instrument developed at Helsinki University of Technology [31]. The phantom uses TiO_{2} particles and an infrared dye to provide scattering and absorption properties similar to that of biological tissue. The homogeneous optical coefficients of the background material were approximately *µ*′_{s}=1 ± 0.1 mm^{-1} and *µ _{a}*=0.01±0.001 mm

^{-1}at a wavelength of 800 nm. The speed of light for the phantom material was

*c*=0.19 mm ps

^{-1}. Two cylindrical inhomogeneities of diameter 9.5 mm and height 9.5 mm were located in the central plane

*z*=0 of the cylinder. The optical properties of the two targets relative to the background were set to (

*µ*,2

_{a}*µ*′

_{s}) and (2

*µ*,

_{a}*µ*′

_{s}), respectively. The geometry of the phantom is shown in Fig. 2. Optical fibres were used to transmit light from the source and to the detectors. In this experiment, 16 source and 16 detector sites arranged in two rings at a spacing of 12 mm were used. Filtering out the measurements from the detectors that are close to a source when the source is used to avoid light saturation, results in a final set of 192 amplitude and phase measurements in our data sets. To suppress amplitude scalings between the experimental and simulated datasets, data was collected in two sets, one with the blobs inside the homogeneous domain g

^{meas}and one where the blobs were absent g

^{ref}. The augmented residual of (17) becomes

Where *A*(**p**
^{ref}) represents the solution of the numerical model without any blobs inside and homogeneous optical properties of *µ*
_{a}=0.0078*mm*
^{-1} and *µ*′_{s}=1.065*mm*
^{-1}

The Finite Element Model used in the reconstruction consisted of an unstructured cylindrical mesh with 22351 nodes and 121500 linear tetrahedral elements. The initial guesses for the unknown blobs were parameterised using Spherical Harmonics of second degree, represented by sets of twelve coefficients. The mapped surface mesh was of 1442 nodes and 2880 linear triangular elements.

The reconstruction results are shown in Fig. 3(e) and Fig. 4(e) for the absorption target and Fig. 3(f) and Fig. 4(f) for the scattering target. Both the absorption and scattering targets were recovered well in terms of general shape as well as localization. In the same figures we show a drawing of the actual positions for the targets as given by the experiment and the solution from the use of a pixel based algorithm as presented in [34]. Since the shape based method constrains any cross talk and artifacts from appearing the reconstructed image produces a much higher contrast image than the pixel based algorithm, which we can see by comparison in Fig. 4 and Fig. 3.

To assess the success of the shape based reconstruction we present in Fig. 5 a set of metrics which we calculated for each iteration of the algorithm. The first one is the normalized residual of the data

which can be seen in Fig. 5(a). The residual is decreased during the first iterations until a local minimum was reached and then with the addition of further degrees of spherical harmonics is allowed to decrease further.

To measure the rate of convergence for the algorithm in the parameter space, we calculated, in Fig. 5(b), the ratio of the volumes of the reconstructed shapes to the volumes of the target shapes for each iteration. To approximate the shape properties of the targets we used cylindrical surface meshes Γ^{target}
_{1} and Γ^{target}
_{2} for the *µ*
_{a} and the *µ*′_{s} targets, based on the assumed sizes and locations. To calculate the volume each object was divided into tetrahedra based on the mapped mesh and the volumes of all the tetrahedra were later added to calculate the total volume of the object. The desired value for this ratio is one and we notice that both the absorption and the scattering objects are performing quite well with the absorption reconstructed object being smaller than the expected target. The main reason for the smaller absorption target volume could be an inaccuracy in the experimental calculation of the “correct” values for both the absorption and scattering for the background and the targets. A study of the relation between the optical parameters and the volume was presented in [25] and were consistent with theoretical results showing an inherent ambiguity between size and contrast for diffusion coefficient [36] and for absorption coefficient [37].

To assess how close the reconstructed shapes were to the expected ones we also calculated the Hausdorff distance between the surfaces of the cylindrical targets and the recovered surface in each iteration in Fig. 5(c) The Hausdorff distance measures the degree of mismatch between the two sets, as the distance of the node *v* of Γ* ^{target}_{p}* that is the farthest away from any node of the evolution surfaces Γ

*at iteration n and is defined as*

^{n}_{p}$H\mathrm{aus}({\Gamma}_{p}^{t\mathrm{arg}\mathrm{et}},{\Gamma}_{p}^{n})=\underset{\mathbf{v}\in {\Gamma}_{p}^{target}}{\mathrm{max}}\underset{\mathbf{v}\prime \in {\Gamma}_{n}^{p}}{\mathrm{max}}\mid \mathbf{v}-\mathbf{v}\prime \mid $

where *p*∈[1, 2]. We notice that even though the image of the reconstructed targets along the *z*=0 plane is accurate the sharp edges on the top and bottom of the cylindrical targets were not reconstructed equally well, as can be seen in Fig. 3(e) and Fig. 3(f). This is an expected result due to the diffusive nature of DOT, together with the positioning of the sources and detectors along the *z*=0 plane that reduces the sensitivity of the the data in the z-direction.

Finally, in Fig. 5(d), the zeroth degree of the spherical harmonics representation, which represents the center of the objects, is used to calculate the distance between the expected center of the targets and the evolution objects in each iteration.

The resulting graph shows that the center of the expected targets was recovered successfully. As a further verification of the proposed algorithm we performed the reconstruction repeatedly starting from three different initial guesses, as can be seen in (media 1), (media 2) and (media 3) from Fig. 6. The shape based method performed quite well recovering the targets in each case as can be seen in the movies linked to this article. This provides evidence that convergence to the desired solution is not heavily depended on the choice of the initial guesses.

## 6. Conclusion

We presented a shape based reconstruction algorithm for diffuse optical tomography. The forward model was based on a Finite Element Method with a novel mesh subdivision algorithm for the mapping from the spherical harmonics coefficients to the absorption and scattering pa-rameter distributions in an unstructured volumetric FEM mesh. We validated the shape based method by reconstructing successfully a scattering and an absorption target from experimental measurements in a cylindrical phantom.

The Finite Element Method proved a good tool for three dimensional shape based reconstructions since it can deal with different topologies without a change in the size of the system matrix and could allow for reconstructions to include piecewise constant regions together with non-homogeneous regions.

The advantage of the method relies on the reduction in the dimension of the inverse problem leading to fast reconstructions, that are free from local voxel artifacts. This method could have many applications in the experimental data reconstructions and could enhance the quality of images in DOT reconstructions.

In this paper we showed that the proposed method works with experimental measurements. There is still work to be done to validate the shape based approach for a clinical application. More specifically, some of the problems that we might face arise from the fact that biological tissue is not, in general, piecewise constant. To initially assess this problem we conducted simulation experiments wherein the background region was not homogeneous but drawn from a Normal distribution with specifed variation. In these tests (results not included in this paper) we found the proposed method tolerated variation up to 10% in the absorption and scatter of the background without significant decrease in the recovery accuracy of the inclusions. Furthermore, the largest parameters contrast in biological tissues is usually expected across the tissue interfaces. The choice of the FEM in the shape-based method was so that the homogeneity of all the regions in the domain is not a requirement and homogenous regions can be combined with non homogeneous in the same problem. This extension is part of the future work towards reconstruction of data from biological tissues. In [25] we have seen that biological tissues can be modeled using a parametric shape based method and in [38] we have seen a combined method that uses a image based method to create an good initial guess for the positions of the different types of tissue and a shape based for the creation of a improved reconstruction.

In addition, since atlases of biological tissue are widely available, and the possibility of multimodality instrumentation becomes more and popular, the inclusion of good initial information for a shape based approach can be achieved.

We believe that the demonstration of the shape-based method able to produce good results for experimental data, is a significant step towards future investigation on the feasibility of the method for clinical applications.

## Acknowledgments

The authors would like to acknowledge D.Sc. Ilkka Nissila for providing the experimental data and Dr. Oliver Dorn for many useful discussions. The work was supported by EPSRC grant EP/E034950/1 and the Academy of Finland (projects 119270 and 213476)

## References and links

**1. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems **15**, 41—93 (1999).
[CrossRef]

**2. **A. G. Yodh and D. A. Boas, “Functional imaging with diffusing light,” Biomedical photonics handbook(CRC Press, 2003).

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

**4. **J. P. Kaipio, V. Kolehmainen, M. Vauhkonen, and E. Somersalo, “Inverse problems with structural prior information,” Inverse Problems **15**, 713–729 (1999).
[CrossRef]

**5. **A. Douiri, M. Schweiger, J. Riley, and S. R. Arridge, “Anisotropic diffusion regularisation methods for diffuse optical tomography using edge prior information,” Meas. Sci. Tech. **18**, 87—95 (2007).
[CrossRef]

**6. **M. Schweiger and S.R. Arridge, “Optical Tomographic Reconstruction in a Complex Head Model Using a priori Region Boundary Information,” Phys. Med. Biol. **44**, 2703–2721 (1999).
[CrossRef] [PubMed]

**7. **C. Panagiotou, S. Somayajula, A. P. Gibson, M. Schweiger, R. M. Leahy, and S. R. Arridge, “Information theoretic regularization in diffuse optical tomography,” J. Opt. Soc. Am. A **26**, no. 5, 1277–1290 (2009).
[CrossRef]

**8. **P. Hiltunen, D. Calvetti, and E. Somersalo, “An adaptive smoothness regularization algorithm for optical tomography,” Opt. Express **16**, 19957–19977 (2008).
[CrossRef] [PubMed]

**9. **S. R. Arridge, O. Dorn, V. Kolehmainen, M. Schweiger, and A. Zacharopoulos, “Parameter and structure reconstruction in optical tomography,” Journal of Physics: Conference Series 125,012001 (2008).

**10. **J. Chang, H.L. Graber, P. C. Koo, R. Aronson, S.-L. S. Barbour, and R. L. Barbour, “Optical Imaging of Anatomical Maps Derived from Magnetic Resonance Images Using Time-Independent Optical Sources,” IEEE Trans. **16**, 68–77,(1997).

**11. **B. W. Pogue and K. D. Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of a priori magnetic resonance imaging structural information,” Opt. Lett. **23**, Is.21, 1716–1718, (1998).
[CrossRef]

**12. **V. Ntziachristos, A. G. Yodh, M. D. Schnall, and B. Chance, “MRI-guided diffuse optical spectroscopy of malignant and benign breast lesions,” **4**, 347–354, (2002).

**13. **M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. **50**, 2837–58, (2005).
[CrossRef] [PubMed]

**14. **X. Intes, C. Maloux, M. Guven, B. Yazici, and B. Chance, “Diffuse optical tomography with physiological and spatial a priori constraints,” Phys. Med. Biol. **49**, N155–63, (2004).
[CrossRef] [PubMed]

**15. **B.A. Brooksby, H. Dehghani, B.W. Pogue, K.D. Paulsen, and KD. “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE Sel. Top. Quantum Electron. **9**, 199–209, (2003).
[CrossRef]

**16. **S. R. Arridge, C. Panagiotou, M. Schweiger, and V. Kolehmainen, “Translational multimodality optical imaging, Chapter 5: Multimodal Diffuse Optical Tomography: Theory,” Artech Press, 101–124 (2008).

**17. **S. R. Arridge and J.C. Schotland, “Optical Tomography: Forward and Inverse Problems”, Inverse Problems **25**, 12, (2009).
[CrossRef]

**18. **O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Problems **22**, R67–R131, (2006).
[CrossRef]

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

**20. **M. Schweiger, S. R. Arridge, O. Dorn, A. Zacharopoulos, and V. Kolehmainen, “Reconstructing absorption and diffusion shape profiles in optical tomography by a level set technique,” Opt. Lett. **21**, no. 4, 471–473 (2006).
[CrossRef]

**21. **N. Naik, R. Beatson, J. Eriksson, and E. van Houten, “An implicit radial basis function based reconstruction approach to electromagnetic shape tomography,” Inverse Problems25, no. 2 (2009).
[CrossRef]

**22. **D. Alvarez, P. Medina, and M. Moscoso, “Fluorescence lifetime imaging from time resolved measurements using a shape-based approach,” Opt. Express **17**, 8843–8855 (2009).
[CrossRef] [PubMed]

**23. **M. Soleimani, W. R. B. Lionheart, and O. Dorn, “Level set reconstruction of conductivity and permittivity from boundary electrical measurements using experimental data,” Inverse Problems in Sc. and Eng. **14**, 193–210 (2006).
[CrossRef]

**24. **M. Soleimani, O. Dorn, and W. R. B. Lionheart, “A narrow-band level set method applied to EIT in brain for cryosurgery monitoring,” IEEE Trans. Biomed. Eng. **53**, 2257–2264 (2006).
[CrossRef] [PubMed]

**25. **A. Zacharopoulos, S.R Arridge, O. Dorn, V. Kolehmainen, and J. Sikora, “Three-dimensional reconstruction of shape and piecewise constant region values for optical tomography using spherical harmonic parametrisation and a boundary element method,” Inverse Problems **22**, 1509–1532 (2006).
[CrossRef]

**26. **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]

**27. **G. Boverman, E.L. Miller, D.H. Brooks, D. Isaacson, F. Qianqian, and D.A. Boas, “Estimation and statistical bounds for three-dimensional polar shapes in diffuse optical tomography,” IEEE Trans. Med. Imaging **27**, 752–765 (2008).
[CrossRef] [PubMed]

**28. **M. E. Kilmer, E. L. Miller, A. Barbaro, and D. Boas, Three-dimensional shape-based imaging of absorption perturbation for diffuse optical tomography,” Appl. Opt. **42**, 3129–3144 (2003).
[CrossRef] [PubMed]

**29. **S.R. Arridge, M. Schweiger, M. Hiraoka, and D.T. Delpy, “A finite element approach for modelling photon transport in tissue,” Med. Phys. **20**, 299–309 (1993).
[CrossRef] [PubMed]

**30. **A. Björck, “Numerical methods for least square problems,” SIAM, (1996).

**31. **I. Nissilä, K. Kotilahti, K. Fallström, and T. Katila., Instrumentation for the accurate measurement of phase and amplitude in optical tomography,” Rev. Sci. Instrum. **73**, 3306–3312 (2002).
[CrossRef]

**32. **T. Moller and B. Trumbore, *Fast, minimum storage ray-triangle intersection*,” J. Graphics Tools **2**, 21–28 (1997).

**33. **R. Kopperman, P. Meyer, and R.G. Wilson, “A Jordan surface theorem for three-dimensional digital spaces,” Discrete Comput. Geom. **6**, 155–161 (1991).
[CrossRef]

**34. **M. Schweiger, S. Arridge, and I. Nissila, “GaussNewton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol. **50**, 2365–2386 (2005).
[CrossRef] [PubMed]

**35. **C.R. Vogel, “Computational methods for inverse problems,” *Frontiers in Applied Mathematics Series*, SIAM, 23 (2002).

**36. **D.J. Cedio-Fengya, S. Moskow, and M.S. Vogelius, “Identification of conductivity imperfections of small diameter by bounday measurements. Continuous dependence and computational reconstruction,” Inverse Problems **14**, 553–595 (1998).
[CrossRef]

**37. **G. Bal, “Optical tomography for small volume absorbing inclusions,” Inverse Problems , **19**, 371–386 (2003).
[CrossRef]

**38. **V. Kolehmainen, S.R. Arridge, M. Vauhkonen, and J.P. Kaipio, “Recovery of constant coefficients in optical diffusion tomography,” Opt. Express **7**, 468–480 (2000).
[CrossRef] [PubMed]