## Abstract

The level set technique is an implicit shape-based image reconstruction method that allows the recovery of the location, size and shape of objects of distinct contrast with well-defined boundaries embedded in a medium of homogeneous or moderately varying background parameters. In the case of diffuse optical tomography, level sets can be employed to simultaneously recover inclusions that differ in their absorption or scattering parameters from the background medium. This paper applies the level set method to the three-dimensional reconstruction of objects from simulated model data and from experimental frequency-domain data of light transmission obtained from a cylindrical phantom with tissue-like parameters. The shape and contrast of two inclusions, differing in absorption and diffusion parameters from the background, respectively, are reconstructed simultaneously. We compare the performance of level set reconstruction with results from an image-based method using a Gauss-Newton iterative approach, and show that the level set technique can improve the detection and localisation of small, high-contrast targets.

© 2010 Optical Society of America

## 1. Introduction

Optical diffusion tomography (ODT) is a non-invasive and portable imaging technique that recovers the optical parameters of a heterogeneous scattering medium from boundary measurements of light transmission in the near-infrared wavelength range. The main applications are in diagnostic medical imaging, where the reconstructed volume distributions of absorption and scattering parameters can be related to physiological states and processes which cannot be imaged by other diagnostic methods [1] such as blood and tissue oxygenation levels, or oxygen uptake and metabolism rates. The applications include monitoring of oxygenation levels in term and preterm infants, brain activation studies [2, 3] and breast tumour screening.

In many biomedical applications the parameter distribution sought by the inverse solver contains sharply defined objects with step-like boundaries. These can for example be the outlines of a tumour or haematoma, or the interfaces between different organs or tissue types, or the boundaries of regions filled with some tracer or marker substance. These interfaces are typically not recovered well in classical image-based reconstruction schemes that recover the parameter values for each voxel, due to the need for regularisation, which may blur sharp edges or flatten the contrast of inclusions. Most regularisation tools penalise variations or gradients in the parameter distribution, yielding over-smoothed reconstructions. However, in many applications the existence of high-contrast inclusions is known a-priori, and it is important to be able to find the shape of inclusions accurately. Whereas the parameter values may vary only slightly within the inclusion and the background medium, often across the interfaces significant jumps occur.

Shape-based reconstruction methods that recover the shape and location of the interfaces provide an alternative to voxel-based reconstruction of the parameter distributions in this case. Explicit shape methods parametrize the interface surfaces between regions, e.g. by a spherical harmonics expansion [4] or an ellipsoidal description [5], and reconstruct the regions by recovering the shape parameters. An alternative approach to shape-based reconstruction is the *level set* technique, where the zero contours *ψ*(**r**) = 0 of a level set function *ψ* defined over the domain Ω define the boundary between regions with a step-like parameter contrast. Level set methods are a class of implicit shape methods, i.e. the topology of the problem is not predefined and can change during the image reconstruction process. An application for a shape-only recovery in electromagnetic tomography has been previously presented in [6]. An example of combined shape and a single physical parameter reconstruction in electrical impedance tomography has been presented in Chung *et al*[7]. Bal and Ren applied the level set approach to an inverse interface problem of recovering singular surfaces of clear regions [8]. Irishina *et al* [9] use the level set method for early detection of breast cancer with microwave imaging. Jacob *et al* used an iterative approach for recovering brain activation images that alternated between a level set step and an image-based step [10]. The techniques for an efficient calculation of gradient directions are given in Tai. *et al* [11]. A recent review is given in Dorn et al. [12]. We have previously discussed the application of a level set method to the simultaneous recovery of absorption and scattering inhomogeneities from simulated model data in a two-dimensional [13] and three-dimensional test problem [14]. In this paper we apply the level set approach to experimental data obtained from phantom measurements, and compare the results with an image-based Gauss-Newton reconstruction scheme.

## 2. Method

Light transport in scattering media can be described by the radiative transfer (RTE) equation [15, 16, 17, 18]. Infrared light transport in most biological tissues is scatter-dominated, where photons contributing to the measurements have undergone multiple scattering events. In this case the RTE can be reduced to the computationally efficient diffusion equation (DE) [19, 20, 21, 22]. Given a scattering domain Ω⊂ℝ_{3}, bounded by ∂Ω, the diffusion equation in the frequency domain is given by

with Robin-type boundary condition at the tissue-air interface ∂Ω:

where *ω* ∈ ℝ^{+} is the modulation frequency, Φ is the photon density field, *c* = *c*
_{0}/*n* is the speed of light in a medium with refractive index *n*, *c*
_{0} is the vacuum speed of light, *q* is a source term formulated as an incident flux, *ζ* is a boundary term which incorporates the refractive index mismatch at the surface, *ν* is the outward normal at ∂Ω, and *κ* and *μ _{a}* are the diffusion and absorption coefficients, respectively. In the following we assume each source

*q*to be a radio-frequency modulated signal at a frequency

_{i}*ω*

_{0}:

where *u _{i}* defines a source intensity profile over ∂Ω centered at source position

**m**

_{i}^{(S)}. The boundary measurements

**M**

*obtained from an experimental data acquisition system are represented in the model by the integral*

_{ij}**y**

*of the outgoing diffuse exitance Γ for source*

_{ij}*i*,

over a detector profile *W _{j}* on the boundary, centered at detector position

**m**

_{j}^{(D)}, where Γ is given by

The model parameters {*x*
_{1}(**r**),*x*
_{2}(**r**)} = {*κ*(**r**),*μ _{a}*(

**r**)} are discretised into a finite-dimensional vector of basis coefficients

**x**= {

**x**

_{1},

**x**

_{2}} ∈ ℝ

^{N}by a suitable basis expansion

with basis functions *b _{k}*, and

**x**

*= {*

_{l}*x*} is the set of coefficients for the expansion of

_{lk}*x*(

_{l}**r**).

In general the basis representations of the volume parameters for the forward and inverse problems may be different. We used a piecewise polynomial basis on an unstructured grid to solve the diffusion equation numerically with a finite element model. For the solution of the inverse model, the parameters were expanded into a tri-linear basis on a regular grid [23]. A linear transformation was applied to map between the two basis representations.

The image reconstruction problem in ODT is generally approached as a model-based optimisation problem, where a model of light transport in scattering media is compared to the measurement data and the parameters of the model are iteratively updated such as to minimise the difference to the measured data. The forward model defined by Equations (1)–(6) can be formally expressed as an operator ℱ,

The reconstruction of the parameter images **x**̂, given a set of measurement data **M**, can then be expressed in terms of the minimisation of an *objective function* 𝓙** _{M}**,

where a common choice of 𝓙 is a regularised least squares cost functional,

where the data misfit is expressed by residual operator 𝓡** _{M}** that describes the difference between model data

**y**and measurement data

**M**

and *Q* is a regularisation functional multiplied by the hyperparameter *τ*.

Standard approaches for solving the optimisation problem Eq. 7 include methods that make use of first order derivative information, such as nonlinear conjugate gradients [24], or second order derivative information, such as Newton-type methods [25, 26, 22, 27]. The reconstruction problem in optical tomography is ill-posed and significantly affected by data noise. The regularisation term *Q* is required to impose additional constraints on the solutions. Typically, *Q* applies local smoothness conditions on the solution. Positivity of the parameter distributions *x _{l}* can be enforced by mapping to logarithmic values.

Where the target image is known to consist of high-contrast inclusions in a relatively flat background, an alternative to conventional pixel-based reconstruction techniques is the utilisation of shape-based inverse methods.

In the shape inverse problem we assume that the domain of interest is divided into several disjoint zones with distinct optical parameter distributions. In the simplest case, we consider Ω to be divided into an *object region S _{l}* and a

*background region*Ω\

*S*for each parameter

_{l}*x*(

_{l}**r**),

*l*∈ {1,2}, where

*S*

_{1}and

*S*

_{2}need not coincide.

If we consider the parameter distribution to be piecewise constant and comprised of two distinct intensity levels, then *x _{l}* is of the form

where *x*
_{l,i} and *x*
_{l,e} are constants. If their values are known, the task of the reconstruction is to recover the shapes *S _{l}*. If

*x*

_{l,i,e}are not or only partially known a-priori, they can be added as unknowns to the reconstruction problem. Typically, the background parameters

*x*

_{l,e}may be assumed to be known, while the inclusion contrasts

*x*

_{l,i}are to be reconstructed.

To define the regions *S _{l}* with the level set technique[6, 28, 29, 30], we introduce sufficiently smooth level set functions

*ψ*such that

_{l}The boundary *D _{l}* = ∂

*S*of object

_{l}*S*is defined by the zero level set of the level set function

_{l}*ψ*. To solve the shape reconstruction problem, we will adopt a time evolution approach [29]. As a consequence,

_{l}*D*and

_{l}*ψ*will be functions of an artificial evolution time

_{l}*t*, i.e.

*D*(

_{l}*t*) = {

**r**:

*ψ*(

_{l}**r**,

*t*) = 0}.

The inverse problem can be stated as follows: For a shape-only reconstruction, find level set functions *ψ _{l}* that minimise the shape least squares cost functional

where we denote 𝓡** _{M}**(

*ψ*) = 𝓡

_{l}**(**

_{M}*x*(

_{l}*ψ*)). For a combined shape and object contrast reconstruction, the object contrast also evolves with time,

_{l}*x*

_{l,i}=

*x*

_{l,i}(

*t*), and in addition to level set functions

*ψ*we find interior parameter values

_{l}*x*

_{l,i}that minimise the least squares cost functional

where we denote 𝓡** _{M}**(

*ψ*,

_{l}*x*

_{l,i}) = 𝓡

**(**

_{M}*x*(

_{l}*ψ*,

_{l}*x*

_{l,i})). We want to derive an evolution equation for the unknown parameters

*ψ*

_{1}(

*t*),

*ψ*

_{2}(

*t*), and optionally

*x*

_{1,i}(

*t*) and

*x*

_{2,i}(

*t*) which reduces (and eventually minimizes) the above defined least squares cost functional. We are therefore looking for forcing terms

*f*

_{1}(

**r**,

*t*),

*f*

_{2}(

**r**,

*t*), and optionally

*g*

_{1}(

*t*), and

*g*

_{2}(

*t*) such that the system

defines a descent flow for (9) or (10), i.e.

Following the derivation given in [13], we can evaluate the derivatives of 𝓙_{M}^{(s)} and 𝓙_{M}^{(s+c)} with respect to the evolution time to be

$$\frac{d{\U0001d4d9}_{\mathbf{M}}^{\left(S+C\right)}}{\mathrm{dt}}=\frac{d{\U0001d4d9}_{\mathbf{M}}^{\left(S\right)}}{\mathrm{dt}}+\sum _{l=1}^{2}\mathrm{Re}\left({g}_{l}{\int}_{\Omega}\left(1-H\left({\psi}_{l}\right)\right){\U0001d4e1\prime}_{l}{\left({x}_{l}\right)}^{*}\U0001d4e1d\mathbf{r}\right),$$

where the symbol Re indicates to take the real part of the complex-valued arguments, 𝓡′* _{l}*(

*x*)

_{l}^{*}is the adjoint Fréchet derivative, and

*δ*(

*ψ*) =

_{l}*H*′(

*ψ*) is the one-dimensional Dirac delta distribution. In the following numerical reconstructions the expressions 𝓡

_{l}*′(*

_{l}*x*)

_{l}^{*}𝓡 are determined using an adjoint scheme [31].

Using the fact that *δ*(*ψ _{l}*) > 0, we can define descent directions

*f*

_{l,d}and

*g*

_{l,d}for

*f*and

_{l}*g*as

_{l}$$\phantom{\rule{1.em}{0ex}}{g}_{l,d}\left(t\right)=-\mathrm{Re}\left({\int}_{{S}_{l}}{\U0001d4e1\prime}_{l}{\left({x}_{l}\right)}^{*}\U0001d4e1d\mathbf{r}\right),$$

for *l* = 1, 2 which satisfy the descent flow condition (12) for 𝓙.

In the presence of data or parameter noise, an additional regularisation scheme must be applied. Data noise may arise from limited signal to noise ratio in the data acquisition system, or systematic noise such as surface coupling losses. Parameter noise in this context is introduced by incomplete knowledge of the background parameter distributions *x*
_{l,e}(**r**) such that *x _{l}* is not in the range of the model ℱ(

*x*).

_{l}Instead of applying an explicit regularisation functional *Q* to the cost term (7), we perform an implicit regularisation by applying a smoothing operator 𝓟 to the forcing terms in Eq. 13 which has the effect of smoothing the boundaries of the shapes during the evolution by projecting the updates towards a smoother subspace[32]. We use a 𝓟 of the form

where Δ denotes the Laplace operator and where the regularization parameters *α* > 0 and *β* > 0 can be chosen freely. Discretizing (11), (13) by a straightforward finite difference time-discretization with time-step Δ*t* > 0 and interpreting *ψ _{l}*

^{(n)}=

*ψ*(

_{l}*t*),

*ψ*

_{l}^{(n+1)}=

*ψ*(

_{l}*t*+ Δ

*t*),

*f*

_{l,d}

^{(n)}=

*f*

_{l,d}(

*t*) and

*g*

_{l,d}

^{(n)}=

*g*

_{l,d}(

*t*), we arrive at the iteration

$${\psi}_{l}^{\left(n+1\right)}={\psi}_{l}^{\left(n\right)}+\U0001d4df\Delta t{f}_{l,d}^{\left(n\right)},$$

$${x}_{l,i}^{\left(n+1\right)}={x}_{l,i}^{\left(n\right)}+\Delta t{g}_{l,d}^{\left(n\right)}.$$

## 3. Results

We present level set reconstructions from data generated by a numerical diffusion light transport model, and from experimental data obtained from a phantom with tissue-like optical parameters. The phantom consisted of a cylindrical domain with homogeneous background and two embedded cylindrical objects, with increased absorption and increased scattering parameters, respectively. To allow a direct comparison between the model and experimental data reconstructions, we simulated the same object geometry, parameter distributions and measurement arrangement in the model data calculations. As a reference, we also performed reconstructions of the model and phantom data using an image-based iterative Gauss-Newton-Krylov scheme [22].

The phantom had a radius 35 mm and height 110 mm, oriented in a reference frame such that Ω= {(*x*,*y*,*z*)|(*x*
^{2} + *y*
^{2})^{1/2} ≤ 35 mm and |*z*| ≤ 55 mm}. The background parameters were homogeneous with *x*
_{1,e}(**r**) = *μ*
_{a,e} = 0.0078mm^{-1} and *x*
_{2,e}(**r**) = *κ _{e}* = 0.31 mm. The two inclusions had cylindrical shape with radius 9.5 mm and height 9.5 mm, located in the central plane of the phantom, as shown in Fig. 1. One of the inclusions was centered at position (-17.5,0,0) and had an increased absorption coefficient of

*x*

_{1,i}= 0.0156mm

^{-1}, the other was centered at position (30.3,-17.5,0) and had a decreased diffusion coefficient of

*x*

_{2,i}= 0.155 mm. The refractive index was homogeneous over the cylinder with

*n*= 1.56.

Sources and detectors were arranged in two rings around the cylinder mantle at altitudes *z* = ±6 mm. Each ring contained 8 sources and 8 detectors, indicated with ‘*x*’ and ‘o’ in Fig. 1. For the simulations, both the source profiles *u _{i}*, and detector profiles

*W*were assumed to be Gaussian with width

_{j}*σ*= 2 mm. The boundary measurements consisted of the logarithmic amplitude and phase shift for a modulation frequency of

*ω*

_{0}= 2

*π*× 10

^{8}Hz.

In addition to measurements **M**
^{(sig)} in the presence of the inclusions, reference data **M**
^{(bkg)} were obtained on a homogeneous object with *x _{l}* =

*x*

_{l,e}. The reference data were then used to generate synthetic data

**M**̃ with

and **M**̃ was used to evaluate the misfit operator (Eq. 8) for the level set reconstructions. Difference reconstructions can eliminate systematic differences between model and experiment. This includes uncertainties in signal amplitude, coupling losses, errors in the outer boundary of the object, or, in general, unknown model errors of any kind. Practical applications often require the reconstruction of parameter differences from measurements taken at different subject states, for example during stimulus and rest, or at different times to follow treatment success, or measurements at different wavelengths to classify the object type.

#### 3.1. Reconstructions from simulated data

Measurements were simulated with a finite element model using a mesh consisting of 104577 nodes and 589824 tetrahedral elements with piecewise linear shape functions. The simulated measurement data were then contaminated with 0.05% additive Gaussian noise.

In addition to the data noise, the background optical parameter distribution was modulated with random Gaussian noise of different magnitude levels, to simulate the presence of unknown background variations. The parameter noise was generated by perturbing each pixel with a normally distributed random value, followed by smoothing over the image. 50 different random realisations of the background variations were generated for each of 3 different noise levels *σ _{p}* = 10%, 20% and 50% of the mean background value. The background distributions were assumed to be identical for signal and reference measurements. Fig. 2 shows cross sections of one sample of absorption and diffusion distribution for each of the noise levels.

**a) Reconstructions for shape only**. In this section, we present results of level set difference reconstructions from simulated data, where object contrast was assumed known and fixed (shape-only reconstruction). Simulated measurement data were generated for the homogeneous background case, as well as for all 50 samples of each of the three background noise levels indicated in Fig. 2. Level set reconstructions were performed for each data set.

Cross sectional images through the reconstructed parameter distributions are shown in Fig. 3. For the cases with random variation of the background distribution, the images show the mean and variance distributions over the 50 samples. It can be seen that the locations and shapes of both objects are recovered very well for all background noise levels, although at higher noise levels the variability of the reconstructed object boundary is slightly increased. Good localisation is also observed in the line profiles through the objects in the reconstructed images shown in Fig. 4, which are identical for all reconstructions up to 20% background noise. Only the 50% parameter noise case shows a small displacement of the objects.

To quantify shape recovery with the level set approach, we define a *shape error ε* as the relative volume of mislabeled regions,

The shape errors for shape-only absorption and diffusion reconstructions are shown with solid lines in Fig. 5 as a function of background noise. It can be seen that the shape errors for the shape-only reconstructions, under the assumption of known object contrasts, are generally smaller than the shape errors for the combined shape and contrast reconstructions, indicating cross-talk between shape and contrast errors. In addition, the shape errors increase with increasing background noise, except for a small minimum at noise level 10% for the diffusion inclusion.

b) **Reconstructions for shape and inclusion contrast** The results of simultaneous reconstructions for shapes and inclusion contrasts *x*
_{l,i} for level set reconstruction from difference data are shown in Fig. 6. Displayed are the cross sections through the image means and variances over the 50 background noise realisations. Contrast evolutions for one sample of each background noise level are shown in Fig. 7, and cross sections through the objects of those samples are shown in Fig. 8. Again, the inclusion locations and shapes are recovered well for all background noise realisations, with comparable results, although a small artefact can be seen in the diffusion image of the reconstruction from homogeneous background. After 400 level set iterations, the reconstruction from homogeneous background data recovers the absorption contrast very well, but slightly underestimates the diffusion contrast. For noisy backgrounds, the recovered contrast values exhibit more variation, although the shape recovery is little affected.

**c) Comparison with Gauss-Newton-Krylov image-based reconstructions** To estimate the benefits of the level set approach for object identification, we compare the results with reconstructions from an image-based reconstruction of parameter values. The solver in this case employs a damped Gauss-Newton approach using a Krylov linear solver (DGN-K) that supports an implicit definition of the Hessian matrix [22]. A total variation (TV) regularisation scheme is chosen. The reconstruction is performed in a 64×64×64 grid of voxels.

Difference parameter reconstructions were performed from data for each of the four background parameter distribution realisations. In all cases, the initial parameter distribution for the DGN-K reconstruction was homogeneous. The results are shown in Fig. 9, and the corresponding cross sections through the objects are plotted in Fig. 10. Similar to the level set difference reconstructions, the DGN-K difference reconstruction results are nearly independent of the background parameter distributions. However, it can be seen that the image-based approach significantly underestimates the contrast for both objects, and is less able to recover the location and shapes of the objects than the level set approach. Due to low contrast recovery, the signal-to-noise ratio in the images is lower than for level sets.

The DGN solver uses an explicit representation of the Jacobian matrix, and therefore requires additional memory storage compared to the level set method. A MATLAB implementation using the TOAST Optical Tomography toolbox required 1.07GB of memory for the level set solution, and 2.74GB for the DGN solution. Similarly, the DGN solver is computationally more expensive per iteration, at approximately 400 seconds compared to the level set solver with 18 seconds per iteration for the shape-only reconstruction, and 42 seconds per iteration for the shape and contrast reconstruction. However, the DGN solver converged after typically 15–20 iterations, and may therefore be faster than the level set solver, which required typically 10^{2} iterations for the shape-only reconstruction, and 4 × 10^{2} iterations for the shape and contrast reconstruction.

#### 3.2. Reconstructions from experimental phantom data

The final set of reconstructions was performed on experimentally acquired measurement data obtained from a phantom with tissue-like optical properties [33]. The phantom used TiO_{2} particles and an infrared dye to provide scattering and absorption properties. Phantom geometry and parameters, as well as the measurement arrangement, corresponded to the simulated data in Section 3.1. The frequency-domain data acquisition instrument used to collect the measurement data was developed at Helsinki University of Technology [34]. Optical fibres were used to transmit light from the source and to the detectors. Phase and amplitude data were used to make the reconstructions [35]. For reconstructions from data differences, a set of reference data from a homogeneous phantom without inclusions was measured. From the measurement data, both level set and DGN-K image-based difference reconstructions were performed.

**a) Level set reconstructions.** Signal and reference data were used to reconstruct parameter differences from data differences. Two reconstructions were performed: (i), shape only, assuming a-priori known background parameters and object contrasts, and (ii) shape and inclusion contrast, assuming known background parameters.

The reconstruction results are shown in the left and middle column of images in Fig. 11. Both reconstructions recover object and shape locations well. The combined shape and contrast reconstruction shows a small crosstalk artefact in the diffusion image located at the position of the absorbing inclusion. It also slightly overestimates the object sizes compared to the shape-only reconstruction, and in turn underestimates the object contrasts, as can be seen in the line profiles of recovered contrast levels in Fig. 12.

**b) Comparison with Gauss-Newton-Krylov image-based reconstructions.** For comparison of the shape-based level set method with a voxel-based method, we performed a reconstruction using the DGN-K solver with the phantom difference data. The right column of images in Figure 11 shows the reconstruction results, displayed with the same gray-scale range as the level set reconstructions. Line profiles through the target objects are shown in Fig. 12. The reconstruction does recover the target locations for the absorption and diffusion inclusions, but underestimates the contrast in particular for absorption. In addition, comparison of the profiles with the level set results shows less well-defined target boundaries and some cross-talk in the diffusion image.

## 4. Conclusions

We have presented a shape-based level set reconstruction method that recovers object shape and position in three-dimensional problem domains. The method simultaneously reconstructs absorption and scattering objects. It can either reconstruct the object geometry, assuming known background parameters and object contrast levels, or it can simultaneously reconstruct shapes and object contrasts.

Level set approaches are particularly suited for localisation problems when the availability of reference data allows to perform a reconstruction of parameter differences. In this case, the effect of background parameter inhomogeneities can be mostly eliminated from the data, and the level set reconstruction can be applied to recover localised parameter perturbations robustly.

Compared to image-based reconstructions of local parameter distributions, the level set approach shows superior performance in the identification and localisation of small, well-defined inclusions of high contrast. Image-based approaches tend to recover such objects with low contrast and diffuse boundaries, in particular where data noise requires the application of a smoothing prior.

We have performed level set reconstructions from experimental phantom frequency-domain difference data. Similar to the case of data generated by the numerical model, included objects were recovered well, and the recovery of object locations and shapes was superior to an image-based reconstruction of the same data.

We have shown that level set shape reconstructions provide a viable alternative to conventional image-based techniques for the reconstruction of high-contrast targets in 3-D imaging problems in diffuse optical tomography.

## Acknowlegments

This work was supported by EPSRC grant EP/E034950/1 and the EC Seventh Framework Programme (FP7/2007-2013) under grant agreement n. 201076.

## References and links

**1. **M. Cope and D. T. Delpy, “System for long term measurement of cerebral blood and tissue oxygenation on newborn infants by near infrared transillumination,” Med. Biol. Eng. Comput. **26**, 289–294 (1988). [CrossRef] [PubMed]

**2. **D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Sig. Proc. Magazine **18**, 57–75 (2001). [CrossRef]

**3. **B. W. Pogue, K. D. Paulsen, C. Abele, and H. Kaufman, “Calibration of near-infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating phantoms,” J. Biomed. Opt. **5**, 185–193 (2000). [CrossRef] [PubMed]

**4. **A. Zacharopoulos, M. Schweiger, V. Kolehmainen, and S. R. Arridge, “3D shape based reconstruction of experimental data in diffuse optical tomography,” Opt. Express **17**, 18940–18956 (2009). [CrossRef]

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

**6. **O. Dorn, E. L. Miller, and C. Rappaport, “A shape reconstruction method for electromagnetic tomography using adjoint fields and level sets,” Inverse Probl. **16**, 1119–1156 (2000). [CrossRef]

**7. **E. T. Chung, T. F. Chan, and X.-C. Tai, “Electrical impedance tomography using level set representation and total variational regularization,” J. Comp. Phys. **205**, 357–372 (2005). [CrossRef]

**8. **G. Bal and K. Ren, “Reconstruction of singular surfaces by shape sensitivity analysis and level set method,” Math. Mod. Methods Appl. Sci. **16**, 1347–1373 (2006). [CrossRef]

**9. **N. Irishina, M. Moscoso, and O. Dorn, “Microwave imaging for early breast cancer detection using a shape-based strategy,” IEEE Trans. Biomed. Eng. **56**, 1143–1153 (2009). [CrossRef] [PubMed]

**10. **M. Jacob, Y. Bresler, V. Toronov, X. Zhang, and A. Webb, “Level-set algorithm for the reconstruction of functional activation in near-infrared spectroscopic imaging,” J. Biomed. Opt. **11**, 064029-1 – 064029-12 (2006). [CrossRef]

**11. **X.-C. Tai and T. F. Chan, “A survey on multiple level set methods with applications for identifying piecewise constant functions,” Int. J. Numer. Anal. Mod. **1**, 25–47 (2004).

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

**13. **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. **31**, 471–473 (2006). [CrossRef] [PubMed]

**14. **M. Schweiger, O. Dorn, and S. R. Arridge, “3-D shape and contrast reconstruction in optical tomography with level sets,” in “First International Congress of the International Association of Inverse Problems (IPIA),” J. Phys.: Conf. Ser. **124**, 012043. Institute of Physics, 2007. [CrossRef]

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

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

**17. **K. Ren, G. S. Abdoulaev, G. Bal, and A. Hielscher, “Algorithm for solving the equation of radiative transfer in the frequency domain,” Opt. Lett. **29**, 578–580 (2004). [CrossRef] [PubMed]

**18. **T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. P. Kaipio, “Hybrid radiative-transfer-diffusion model for optical tomography,” Appl. Opt. **44**, 876–886 (2005). [CrossRef] [PubMed]

**19. **K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys. **22**, 691–701 (1995). [CrossRef] [PubMed]

**20. **J. C. Schotland, “Continuous-wave diffusion imaging,” J. Opt. Soc. Am. A **14**, 275–279 (1997). [CrossRef]

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

**22. **M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol. **50**, 2365–2386 (2005). [CrossRef] [PubMed]

**23. **M. Schweiger and S. R. Arridge, “Image reconstruction in optical tomography using local basis functions,” J. Electron. Imaging **12**, 583–593 (2003). [CrossRef]

**24. **S. R. Arridge and M. Schweiger, “A gradient-based optimisation scheme for optical tomography,” Opt. Express **2**, 213–226 (1998). [CrossRef] [PubMed]

**25. **R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation,” Opt. Express **4**, 353–371 (1999). [CrossRef] [PubMed]

**26. **A. D. Klose and A. H. Hielscher, “Quasi-Newton methods in optical tomographic image reconstruction,” Inverse Probl. **19**, 387–409 (2003). [CrossRef]

**27. **J. Nocedal and S. Wright, *Numerical Optimization*, Springer Series in Operations Research and Financial Engineering (Springer, 2006), 2nd ed.

**28. **S. Osher and R. Fedkiw, *Level Sets and Dynamic Implicit Surfaces* (Springer, New York, 2003).

**29. **F. Santosa, “A level-set approach for inverse problems involving obstacles,,” in “ESAIM: Control, Optimization and Calculus of Variations,”, vol. 1 (1996), vol. 1, pp. 17–22.

**30. **J. A. Sethian, *Level Set Methods and Fast Marching Methods* (Cambridge University Press, 1999).

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

**32. **P. González-Rodriguez, M. Kindelan, M. Moscoso, and O. Dorn, “History matching problem in reservoir engineering using the propagation back-projection method,” Inverse Probl. **21**, 565–590 (2005). [CrossRef]

**33. **M. Firbank and D. T. Delpy, “A design for a stable and reproducible phantom for use in near infrared imaging and spectroscopy,” Phys. Med. Biol. **38**, 847–853 (1993). [CrossRef]

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

**35. **I. Nissilä, T. Noponen, K. Kotilahti, T. Tarvainen, M. Schweiger, L. Lipiäinen, S. R. Arridge, and T. Katila, “Instrumentation and calibration methods for the multichannel measurement of phase and amplitude in optical tomography,” Rev. Sci. Instrum. **76** (2005). Art. no. 044302. [CrossRef]