1 Introduction

In an optical diffusion tomography experiment, light from a near infrared laser source is guided to the surface ∂Ω of a highly scattering object Ω by optic fibers, and the transmitted light is measured on ∂Ω by using detecting optic fibres and light sensitive detectors. The objective is to estimate the diffusion and absorbtion coefficients (κ, µ _a) within Ω based on the set of transmission data. The reconstuction of (κ, µ _a) is a nonlinear ill-posed inverse problem.

In this study we consider the reconstruction of κ and µ _a in the case of an object domain Ω⊂ℝ² which is known to consist of L+1 disjoint tissue regions A_k

see Fig. 1. The regions A_k are assumed to be bounded by smooth closed boundary curves {C_ℓ , ℓ=1,…, L} and have constant diffusion and absorbtion coefficients {κ_k , µ _a,k, k=0,1,…, L} within the regions. The outer boundary of the background region A ₀ is ∂Ω.

 

Fig. 1. An example of a piecewise constant object domain Ω=∪ _k , A_k with constant coefficients {κ_k , µ _a,k, k=0, 1, 2, 3}. The outer boundaries of regions A ₁, A ₂ and A ₃ are denoted by C ₁, C ₂ and C ₃, respectively. The background region A ₀ has ∂Ω as the outer boundary.

Download Full Size | PDF

We have recently proposed [1, 2] an approach for the reconstruction of the shapes and locations of the region boundaries {C_ℓ } and the values {κ_k , µ_ak } within the regions {A_k } based on the transmission data on ∂Ω. In the proposed method, the shapes and locations of the boundaries {C_ℓ } are approximated with a set of shape coefficients γ. The FEM discretisation of the forward model, which is the diffusion approximation (DA) to the radiative transfer equation (RTE), was formulated as a mapping from the set of parameters γ, κ=(κ ₀,…, κ_L )^T and µ _a=(µ _a,0,…, µ _a,L)T to the data on the boundary ∂Ω. The inverse problem was then stated as an optimisation problem in which we sought to minimise the residual norm between the measured and computed data.

The performance of the method was evaluated using noisy simulated frequency domain data. In [2] we assumed that the number of regions was known a priori. The method was shown to work well when the initial estimate for the parameters was relatively good in the sense that the location and size of the different tissue regions {A_k } were approximately known a priori. If the location of the regions in the initial estimate were severely erroneous the method could fail to find some regions.

In this paper we show how the shape estimation method can be initialised in the absense of a priori information about the number and locations of the tissue regions. We propose a three stage approach to the recovery of piecewise constant coefficients. As the first stage, constant coefficients ($\tilde{\kappa }$ ₀, $\tilde{\mu }$ _a,0) are fitted to the whole domain Ω using nonlinear least squares. This estimate is used as the initial estimate for the second stage, in which conventional pixelwise images are reconstructed.

In the second stage the inverse problem is stated as Tikhonov regularized problem. As the regularizing penalty functional we use an approximation for the total variation norm of the solution [3, 4, 5]. The regularizing functional is also called a prior due to its identification with the prior distribution in the statistical interpretation of inverse problems [4]. From the statistical point of view, the total variation norm is a feasible choice as a regularizing penalty functional when the images are assumed to consist of a few different levels of {κ_k , µ _a,k} with relatively short boundary lines {C_ℓ }. This is because the total variation penalty functional tends to favor such coefficient distributions which have domains with short boundaries and almost constant values within each region. Once approximate convergence of the second stage is achieved, the initial configuration of the boundaries {C_ℓ } can be found from the pixelwise reconstructions using standard image processing techniques such as edge detection. An example is given in Section 4.

In the third stage, the final images are reconstructed using the shape estimation method that was proposed in [1, 2]. To improve the convergence speed of the the earlier approach, we now also use line search. The performance of this three stage inversion scheme is evaluated by simulations in which the optical parameters (κ, µ _a) are relevant for medical applications of optical tomography.

The rest of this paper is organised as follows. In Section 2 we discuss briefly the parametrisation γ of the region boundaries {C_ℓ } and the diffusion approximation to the RTE. In Section 3 a short stage-by-stage description of the new three stage inversion approach is given and in Section 4 the performance of this approach is evaluated using noisy simulated data. In Section 5 we give conclusions and some suggestions for future work.

2 Forward Model

2.1 Representation of the boundaries {Cℓ}

In this study the domains {A_k } ∊ Ω are assumed to be simply connected and disjoint (i.e. C_i ∩C_j =∅, i≠j), see Fig. 1. We also assume that the outer boundary of the region A ₀, that is, ∂Ω is known and has the property ∂Ω∩C_ℓ =∅ ∀ℓ. If the outer boundaries {C_ℓ , ℓ=1, 2,…,L} of the regions {A_ℓ } are sufficiently smooth, they can be approximated in the form

where θ_n are periodic and differentiable basis functions with period 1. In this study we express both coordinates of the boundary curves as Fourier series with respect to the curve parameter s, that is, we use basis functions of the form

where ϕ^s is the phase of the curve parameter and α denotes either x or y. Furthermore, let γ denote the vector of the shape coefficients, that is,

For more details on the properties of the boundary representation (2) and (3), see [1, 2].

2.2 Diffusion approximation to the RTE

Consider the following experiment: S optic fibers are placed on the source positions ε_j ⊂∂Ω on the boundary of the highly scattering (µ′ _s≫µ _a) object Ω and M optic fibers are placed in the detector positions ζ _i ⊂∂Ω. Light from an intensity modulated near infrared laser source is guided to the object via one of the source fibers at ε_j , and the amount of transmitted light is measured on all the detector locations ζ _i ,i=1,…,M. Then this process is repeated for all S source locations.

A well established model for the above experiment is the diffusion approximation (DA) of the radiative transfer equation [6, 7]. In the case of an intensity modulated light source, the diffusion approximation can be expressed as

where Φ _j is the photon density µ _a is the absorption coefficient [mm^-1], κ=(3(µ _a+µ′ _s))^-1 [mm] is the diffusion coefficient, µ′ _s is the reduced scattering coefficient [mm^-1], c is the speed of light in the medium, q0 _j is the distribution of internal light sources, g _{S j} is the distribution of boundary sources, v is the boundary normal direction, ν is a coefficient due to the mismatch of the refractive indices in Ω and the surrounding medium and the subindex j denotes that the object is illuminated from the source location ε_j .

Two widely used source models within the DA framework are the collimated source (CS) model and the diffuse boundary source (DS) model [8]. In the case of the CS model g _Sj becomes zero, and in the case of the DS model qo _j becomes zero.

The relation between the photon density and the measured exitance at the measurement location ζ _i ⊂∂Ω can be expressed as

For further details on the diffusion approximation and source models, see e.g. [6, 9, 10, 7, 8, 11].

In this study, the discretisation of the forward model (5–6) is based on the finite element method. The details of the discretisation are given in the references [1, 2, 12, 8]. In the following discussion, the vector of image parameters will be denoted by f. The FEM based map P : f↦>z which maps the image parameter vector f to the data vector z on ∂Ω, will be denoted by

The form and size of the parameter vector f in each step of our three stage inversion approach will be given in the next section.

In this paper the data vector is of the form

where Γ _j (ω)=(Γ_1,j(ω),…, Γ_M,j(ω)) is the vector of the complex measurements for the j:th source [13]. This splitting of the data vector means that the discrete mapping (7) ensures that the adjoint operator gives real updates to the solution.

3 Inverse Problem

3.1 Reconstruction of the best fitting constants $\tilde{\mu }$_a,0 and $\tilde{\kappa }$₀

In the first stage of our inversion approach we search for constants ($\tilde{\kappa }$ ₀, $\tilde{\mu }$ _a,0) which minimise the output least squares (LS) error. The dynamic range of the data is very large, and furthermore the real and imaginary parts have quite different scales. We therefore use a L ₂ norm weighted by the inverse of the data diag($z_{\text{meas}}^{-1}$). This type of norm is sometimes called the Rytov Approximation in scattering theory where it is interpreted as a scattered wave normalised by an incoming wave. From a statistical point of view it is equivalent to assuming a diagonal data covariance matrix with the variance of each measurement proportional to the square of the signal and with no correlation between different measurement positions or even between the real and imaginary parts of a single measurement. Although the statistical justification is questionable, since Poisson statistics are arquably more appropriate [14], from an optimisation theory point of view it can be thought of as simply a pre-conditioner for the inverse problem to improve numerical stability [15]. Given this approach, the problem can be stated as

where the parameter vector is f=($\tilde{\kappa }$ ₀, $\tilde{\mu }$ _a,0)^T ∊ ℝ×ℝ. The inverse problem (9) is usually stable (well-posed) due to the small number of unknowns and thus the problem can be solved without additional regularization.

A well-established method for the solution of stable nonlinear optimization problems is the Gauss-Newton method, which can be written formally as

where J̃_(l)=diag($z_{\text{meas}}^{-1}$)J_(l) is the (scaled) Jacobian matrix of P(f ^(l)) and z̃ ^(l)=$z_{\text{meas}}^{-1}$(z _meas-P(f ^(l)).

The motivation behind this stage is to find a good initial estimate for the second step which is the computation of conventional pixelwise reconstructions with the total variation penalty functional.

3.2 Reconstruction in a local pixel basis with the total variation regularization

Once the best matching constants are found, we search an approximative configuration of the regions {A_k } by carrying out reconstruction in a local pixel basis. At this stage, the domain Ω is divided into N_p disjoint pixels Ω _m . In the pixel basis, the functions (κ, µ _a) are approximated in the form

where x_m is the characteristic function of pixel Ω _m . Within the approximation (11), the functions (κ, µ _a) are identified with the vectors ${\mathit{\mu }}_a={({\mu }_{a,1},\dots ,{\mu }_{a,N_p})}^T∊{ℝ}^{N_p}$, $\mathit{\kappa }={({\kappa }_1,\dots ,{\kappa }_{N_p})}^T∊{ℝ}^{N_p}$ and the parameter vector for the inverse problem becomes

At this stage, the ill-posed inverse problem is regularised by using the generalised Tikhonov regularisation. The regularised least squares problem becomes

where ψ(f)>0 is regularizing penalty functional. In this study, ψ(f) is of the form

where TV(·) refers to the total variation norm [3, 4]. Here ${\beta }_{\kappa },{\beta }_{{\mu }_a}$ are constant hyper-parameters whose choice is disscussed in Section 4. With the piecewise constant basis for the coefficients (11), the TV norm can be expressed for both κ,µ _a as

where α is either κ or µ _a, d_j is the length of the jth edge between the adjacent pixels ${\Omega }_{i_1^j}$ and ${\Omega }_{i_2^j}$ and ${\mathbf{\Delta }}_j∊{ℝ}^{N_p}$ is the vector

${\mathbf{\Delta }}_j={(0,\dots ,\stackrel{\left(i_1^{\mathit{j}}\right)}{1,}0,\dots ,0,\stackrel{\left(i_2^{\mathit{j}}\right)}{-1,}0,\dots ,0)}^T$

The choice of the regularizing penalty functional (14) is based on the assumed properties of the true images (κ,µ _a): It is well known that the total variation norm is a feasible prior for piecewise constant images which consist of a few constant levels {κ_k ,µ _a,k} with relatively short well defined boundary lines {C_ℓ } [3, 4]. While the TV prior has found much attention in impedance tomography, it seems to have been used in optical tomography previously only in [5], where positive correlation of the prior distributions of (κ, µ _a) was implicitly assumed, and only the special case of identical contrast in absorption and scattering was examined.

The problem (13) is solved using the Gauss-Newton method. Here, we face the difficulty that the total variation penalty functional is non-differentiable. To overcome this problem, we approximate the L ₁ norm in (15) with the smooth function

where τ>0 is a fixed small valued parameter adjusting the accuracy of the approximation [4]. Another additional feature in the algorithm is due to the fact that (κ,µ _a) are strictly positive functions. To avoid negative values, we use an interior point search and augment the functional in (13) with a penalty term of the form

where {ϛl} is a sequence of decreasing positive coefficients (i.e. ϛl→0 as l→∞), $f_m^{\left(l\right)}$ are the elements of the vector f ^(l) and l denotes the iteration index [16]. With these modifications, the Gauss-Newton algorithm becomes

where ${\mathit{g}}_{{\mathrm{TV}}_{\tau }}^{\left(l\right)}$ and $H_{{\mathrm{TV}}_{\tau }}^{\left(l\right)}$ are gradient and Hessian of the approximated total variation functional ($g_W^{\left(l\right)}$ and ${\mathrm{H}}_W^{\left(l\right)}$ respectively for the functional W(f ^(l)). The details for the computation of the map P(f ^(l)) and the Jacobian matrix J_(l) in equations (10) and (18) can be found in [7]. Details for the computation of ${\mathit{g}}_{{\mathrm{TV}}_{\tau }}^{\left(l\right)}$ and $H_{{\mathrm{TV}}_{\tau }}^{\left(l\right)}$ can be found in [17, 18]. It should be noted that the effect of the functional (17) vanishes as l→∞ and the algorithm (18) converges to the inequality constrained (i.e. f>0) minimum of the functional (13) [16].

The initial estimate for the final stage, which is the actual shape estimation method, can be extracted from the local pixel base reconstructions by using a number of standard image processing tools. An example of feasible image processing steps is given in Section 4.

3.3 Shape estimation

At the final stage of our approach, we utilise a modified version of the shape estimation method which was proposed in [2]. The goal in the method is to find the representation of the boundary shapes γ, and the values κ=(κ ₀,…,κ_L )^T, µ _a=(µ _a,0,…µ _a,L)^T which give the best match to the data in the output LS sense. Thus, the problem can be stated as

where the parameters are

If the number of parameters is not unreasonably high, the inverse problem is not illposed and regularization is not needed. However, to secure numerical stability during the iteration we use a Levenberg-Marquardt type method, which can be formally written as

where d ^(l) is the update direction, s ^(l) is a step length parameter and λ is a stabilisation parameter which is used for the regularization of the search direction d ^(l) It should be noted that in spite of the superficial resemblance of (22) and ordinary Tikhonov regularization, the asymptotic estimate f ^(∞) does not depend on λ - apart from a possibility of getting stuck to a potential local minima. Once the update direction d ^(l) has been computed, the optimal step length parameter s ^(l) can be found by any feasible line search algorithm. We use here the following simple algorithm. Let

denote the one-dimensional output error functional for the line search. To find the minimiser of Ξ_(l)(s), we first search adaptively a set of three points {s ₁, s ₂, s ₃} which are concentrated around the minimum of Ξ_(l)(s) such that 0≤s ₁<s ₂<s ₃ and Ξ_(l)(s ₂)≤min{Ξ_(l)(s ₁), Ξ_(l)(s ₃)}. Then a second order polynomial is fitted to the data {s_i ,Ξ_(l)(s_i )} and the minimum point of the fitted polynomial is used as an approximation for the optimal step length s ^(l)=arg min _s Ξ_(l)(s).

In [2], the optimisation problem (19) was solved in rescaled parameter space in order to avoid numerical instabilities which are due to the different magnitudes of the image parameters γ, κ, µ _a. In this study, we use a different parameter scaling scheme: Instead of scaling the parameters by the initial estimate f ⁽⁰⁾, the parameters are rescaled by $\tilde{\gamma }$ =γ/max(|γ ⁽⁰⁾|), $\tilde{\kappa }$=κ/max(|κ ⁽⁰⁾|) and $\tilde{\mu }$ _a=µ _a/max(|${\mathit{\mu }}_{\mathrm{a}}^{\left(0\right)}$), and the optimisation is performed in the scaled parameter space. This choice yields a numerically more stable problem and speeds up the convergence noticeably.

The details for the computation of the map P(f) and the Jacobian J in (22) can be found in [1, 2]. We note here only that the phase ϕ^s of the curve parameter s in equation (3) has to be fixed when applying equations (21)–(22). If the phase (ϕ^s is not fixed we obtain a one-dimensional null space N_ℓ with respect to each curve C_ℓ in the block J_γ of the Jacobian J.

4 Results

The synthetic measurement system that we use in this paper consists of 16 sources and 16 detectors placed in equiangular positions with the order ε ₁, ζ₁, ε₂, ζ₂,…, ε ₁₆, ζ₁₆ on the boundary ∂Ω of a circular domain Ω which has the radius of 25mm. The optical coefficients within Ω were chosen in the range of interest in medical imaging [7]. For the absorption coefficient µ _a we used values in the range 0.025mm^-1→0.05mm^-1 and for the reduced scattering coefficient µ′ _s we used values in the range 2mm^-1→ 4mm^-1. The diffusion coefficient is computed as κ=(3(µ _a+µ′ _s))^-1. The tissue regions {A_k } were bounded by trigonometric boundaries.

The true target images are shown in the top row of Fig. 2. The regions A ₁ and A ₂ have significant contrast in both the diffusion and absorption coefficients with respect to the background A ₀. The region A ₃, which is located on the lower right side of Ω, has significant contrast in the absorption coefficient and only a weak contrast in the diffusion coefficient with respect to the surrounding region A ₀.

In the generation of the synthetic data with the FEM model we discretised the object domain Ω to P=8600 disjoint triangular elements with D=4429 vertex nodes. The simulated data was computed using a single frequency ω=300MHz and a multiplicative random noise with standard deviation of 1% of the measured signal was added to the simulated data [14]. In the image reconstruction process we used a FEM model in which Ω was divided to P=5776 triangular elements with D=3017 vertex nodes.

In the recovery of the best matching constants we used the initial guess ($\tilde{\kappa }$ ₀, $\tilde{\mu }$ _a,0)=(0.1104,0.0200). The convergence of the iteration was verified by monitoring the output residual norm and the behaviour of the image parameters. The iteration (10) reached full convergence in 12 steps, having minimum at ($\tilde{\kappa }$ ₀, $\tilde{\mu }$ _a,0)=(0.1642,0.0301).

The second stage, which is the reconstruction into the local pixel basis, was started from the constant values (κ ₀, µ _a,0)=(0.1642,0.0301). The domain was discretised into 13 concentric layers of generic quadrilateral pixels, the total number of pixels being N_p =496. Each of these pixels is a union of several triangular FE elements. This disretisation is essentially identical to the “Joshua tree” in [19].

In general, there is not any automatic methods for the choice of two regularization parameters β_κ and ${\beta }_{{\mu }_a}$ in equations (13)–(14), and thus these have to be chosen manually. In the proposed approach the parameters β_κ and ${\beta }_{{\mu }_a}$ are chosen large enough to produce images with flat details and adequately well-defined boundaries. This is adequate since the purpose of this stage is only to produce images which serve as a basis for the initialisation of the boundaries. The initial values of the coefficients are not as crucial as feasible approximate boundaries. We note that the feasible values of β_κ and ${\beta }_{{\mu }_a}$ depend on imaging parameters such as geometry of the domain Ω, data types, noise levels of the instrumentation and the number of pixels in the mesh. For example, for the given imaging parameters we have used fixed values β_κ =50 and ${\beta }_{{\mu }_a}=150$.

The second stage includes the fixed hyperparameter τ in equation (16) and the sequence {ϛl} of the positivity constraint parameters. Consider first the choice of the fixed parameter τ. The function h_τ (t) has the property that h_τ (t)→|t| as τ increases. However, as τ increases the second derivatives of h_τ (t), which are needed in the Hessian matrix ${{\rm H}}_{{\mathrm{TV}}_{\tau }}^{\left(l\right)}$, tend to become very large, leading to numerical instabilities in the iteration (18). Thus, a tradeoff between the accuracy of the approximation (16) and the numerical stability of the iteration (18) has to be made. We have used the value τ=2.5. The choice of the decreasing sequence {ϛl} in iteration (18) is not crucial as long as the positivity constraint remains intact since the asymptotic estimate is not a affected by the positivity constraint [16]. For the given geometry, we have used the sequence {_ϛl+1=0.5 _ϛ,l l≥1} with the initialisation _ϛ0=10^-4. The essential convergence of the iteration (18) was reached in 5 steps. The resulting images are shown in the second row of Fig. 2.

 

Fig. 2. A test case with three boundaries {C_ℓ }. Region A ₁ is located up, A ₂ is located on the lower left side of Ω and A ₃ on the lower right side. Left column shows images of κ([mm]) and right column shows images of µ _a ([mm^-1]). Rows from top to bottom: True distributions, pixelwise reconstructions with the approximated TV prior, initial estimates for the shape estimation and the bottom row shows the final estimates with the shape estimation method. The animation (349kB) shows the progress of the inversion approach step by step.

Download Full Size | PDF

The initial estimate for the shape estimation method was extracted from the pixelwise reconstructions by standard image processing techniques. To make the image processing possible by standard tools, the pixelwise reconstructions were mapped from the “Joshua tree” discretisation to a regular N_r ×N_r grid in the domain G=[-25mm, 25mm]×[-25mm, 25mm] ⊂ℝ². The number N_r was chosen such that the resolution of the regular grid was approximately equal to the resolution of the “Joshua tree” mesh. The left hand image in Fig. 3 shows the reconstructed µ _a in the rectangular mesh. For the grid points which lay outside Ω, we attached values from a pixel on the outermost layer of the Joshua tree mesh. Next, the mapped images of (κ, µ _a) were filtered with a band pass filter and then the edges were detected by looking for zero crossings from the filtered images [20, 21]. The edge detection program returned binary images, having ones at the points where edges were detected, and zeros elsewhere. The right image in Fig. 3 shows the edges that were found from the µ _a reconstruction in the left image of Fig. 3. Based on the reconstructed images and the detected edges, the distinct regions {A_k , k=1, 2, 3} in addition to the background A ₀ were clearly detectable. These regions were marked and boundary curves were fitted to the detected boundary pixels by the LS method, see the animation. These boundary curves and the average values of the coefficients within these regions were then used as the initial guess for the final stage. The initial guess for the final stage is shown in third row of Fig. 2. The initial values of the coefficients {κ_k , µ _a,k, k=0,…,3} can be found in Table 1.

 

Fig. 3. Left: The reconstructed µ _a in a regular grid. Right: The edges that were detected from the µ _a image.

Download Full Size | PDF

Consider the third stage of the proposed approach. In the textbook Levenberg-Marquardt method [22] an iterative search for optimal λ is used at each iteration of the algorithm (21)–(22). However, in the proposed approach we use a fixed λ which is chosen manually such that the algorithm converges. The convergence of the iteration (21)–(22) is verified by monitoring the output error norm and the values of the parameters f ^(l). We note that the choice of λ is not very crucial since the asymptotic estimate

Table 1. Values of {κ_k , µ_a,k} in the initial and in the final estimates of the shape estimation in Fig. 2. The initial estimate is shown in the 3rd row in Fig. 2, and the final estimate in the bottom row, respectively. The dimension of κ_k is [mm] and the dimension of µ_a,k is [mm^-1].

View Table

does not depend on λ. The feasible range for λ depends on the norm of the Jacobian and thus also on the imaging parameters. For the given imaging parameters, we have used typically the value λ=1. Starting from the initial estimate in the third row of Fig. 2, the iteration (21)–(22) converged in 26 iterations. The final estimates for (κ, µ _a) in the case of this example are shown in the bottom row of Fig. 2. The error statistics of the initial and reconstructed values of {κ,_k , µ _a,k} in the third stage are given in Table 1.

As can be seen from Fig. 2, the shapes of the boundaries of the regions A ₁–A ₃ are found with relatively good accuracy. Also the diffusion and absorption coefficients {κ_k , µ _a,k} are found with good accuracy resulting in final parameter estimation error norms in the range ~0–4%, see Table 1.

5 Conclusions

Based on the shape estimation method, which was proposed in [2, 1], we proposed a new three stage scheme for the reconstruction of piecewise constant domains in optical diffusion tomography. In the proposed approach, the shape estimation method is initialised by a local pixel basis reconstruction that is carried out by generalized Tikhonov regularization with a total variation type prior. Standard image processing methods are then used to obtain initial approximations for the boundaries and the coefficients.

The performance of the new approach was evaluated using noisy simulated frequency domain data. It was previously shown that the actual boundary recovery approach is able to yield good estimates for both the boundaries and the coefficients once the initial guess for the iteration was good enough. It was now shown that the proposed approach for the initialization yield the initial configuration with adequate accuracy. Although the performance of the method was evaluated with frequency domain data, the method can be extended to time resolved data types in a straightforward manner. The performance of the method using different data types is a topic of further research.

Although the method assumes that the domains are piecewise constant, it may prove to be useful also in such cases in which the the diffusion and absorption may have small spatial variations within the tissue regions {A_k } and the largest contrast exist over the boundaries {C_ℓ }. The latter approach has also been discussed in [23] wherein an alternative shape estimation method has been investigated. In that case, B-splines with a finite number of control points were used rather than trigonometric functions, and the space of control points was searched with a local greedy strategy and explicit perturbation of a linear forward model.

The method in the form that it was presented here is applicable to 2D situations. In this study we assumed a circular domain, although, in principle, it can also be adopted to other geometries, such as the 2D reflection geometry investigated in [23], by minor modifications to the FEM based forward operator P(f). In that case it may be necessary to include shape regularisation methods as required in [23] although such methods were not needed here.

However, it is well known that optical diffusion tomography should be considered a 3D problem [24]. The validity of a 3D FEM model for optical diffusion tomography has been shown in simulations in [24] and using measured data in [25].

In principle the method presented here can be extended to 3D. The total variation prior has previously been considered in the case of 3D impedance tomography in [18]. The extension of this computational scheme to optical tomography is straightforward as are the image processing operations. In order to utilize the shape recovery we need a shape representation for surfaces, and characterize the inverse problem as the recovery of parameters of these surfaces. A natural choice is spherical harmonics, which have been successfully used for the representation of anatomical structures in other areas of medical imaging [26, 27]. This extension is the obvious next step in the boundary recovery problem.