3D shape based reconstruction of experimental data in Diffuse Optical Tomography

Athanasios Zacharopoulos; Martin Schweiger; Ville Kolehmainen; Simon Arridge

doi:10.1364/OE.17.018940

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 :

- \nabla D (r) \cdot \nabla Φ (r, ω) + μ_{a} (r) Φ (r, ω) + \frac{i ω}{c} Φ (r, ω) = q (r, ω)

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

D = \frac{1}{3 (μ_{a} + μ_{s}^{'})},

Where µ′_s denotes the reduced scattering coefficient. The appropriate boundary condition for our problem is of the Robin type

Φ (r, ω) + 2 ζ D (r) \frac{\partial Φ (r, ω)}{\partial ν} = 0, r \in \partial Ω

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:

y (r) = M Φ = - D (r) \frac{\partial Φ (r; ω)}{\partial ν}, r \in \partial Ω

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

A : 𝒫 \to 𝓩, A (p) = y, p \in 𝒫 and y \in 𝓩

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_i. 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 T_i nodes and a set of associated basis functions υ_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 Φ={ϕ_i} is the h-dimensional vector of basis coefficients and υi are linear basis functions with limited support over the elements that have the node T_i 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.

(K (D) + C (μ_{a}) + \frac{1}{2} α E + i ω B) Φ = q

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

K_{i j} = \int_{Ω} D (r) \nabla υ_{i} (r) \cdot \nabla υ_{j} (r) d r

Equation (6) is formally solved by matrix inversion:

Φ = {(K (D) + C (μ_{a}) + \frac{1}{2} α E + i ω B)}^{- 1} q

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 Ω₀=Ω\∪^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

r ∣ Γ_{ℓ} = {\begin{array}{c} x_{Γ_{ℓ}} (ϑ, φ) = Σ_{l = 0}^{W} Σ_{m = - l}^{l} C_{l, x}^{ℓ, m} {\tilde{Y}}_{l}^{m} (ϑ, φ), \\ y_{Γ_{ℓ}} (ϑ, φ) = Σ_{l = 0}^{W} Σ_{m = - l}^{l} C_{l, y}^{ℓ, m} {\tilde{Y}}_{l}^{m} (ϑ, φ), \\ z_{Γ_{ℓ}} (ϑ, φ) = Σ_{l = 0}^{W} Σ_{m = - l}^{l} C_{l, z}^{ℓ, m} {\tilde{Y}}_{l}^{m} (ϑ, φ), \end{array}

with expansion coefficients {C^m_l}. Here, the basis functions Ỹ^m_l(ϑ,φ) are defined as

{\tilde{Y}}_{l}^{m} (ϑ, φ) ≔ {\begin{array}{c} R e [Y_{l}^{m}] (ϑ, φ), when m \leq 0, \\ I m [Y_{l}^{m}] (ϑ, φ), when m > 0, \end{array}

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

γ = {γ_{k}^{ℓ}} = {C_{l, x}^{ℓ, m}, C_{l, y}^{ℓ, m}, C_{l, z}^{ℓ, m}}, with l = 1, \dots W, m = - l, \dots l

such that the new index k ranges over k=1,⋯3W ². {γ^ℓ_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

μ_{a} = \sum_{ℓ = 0}^{L} μ_{a, ℓ} χ_{ℓ} (r), D = \sum_{ℓ = 0}^{L} D_{ℓ} χ_{ℓ} (r)

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:

K_{ij} = \sum_{ℓ = 0}^{L} \int_{supp (υ_{i} υ_{j}) \cap Ω_{ℓ}} D_{ℓ} (r) \nabla υ_{i} (r) \cdot \nabla υ_{j} (r) d 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_i and T_j nodes.

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_i,y_i,z_i), of the finite element discretisation are inside a given closed surface Γ_ℓ(θ,ϕ). 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_i inside the bounding box with the surface. We discetise the surface Γ_ℓ(θ,ϕ) 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^ℓ_k and nodes S^ℓ_m. We draw the segment R_i={(x_i,y_i,z_i), (x_i,y_i, z_min)}, where z_min is the lower z-coordinate of the volume and check for intersections with each of the triangles G^ℓ_k in the surface approximation Γ_ℓ. 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_i along the xy-plane lies inside. For the triangles satisfying this check we further check for intersections with the line segment R_i, using the [32] ray-triangle intersection algorithm.

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_n of the mesh using their respective nodes. We mark an element as 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 Γ_ℓ.

I (Γ_{ℓ}) = {H_{n} ∣ H_{n} \cap Ω_{ℓ} \neq, and H_{n} \cap Ω'_{ℓ} \neq Ø}

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 ₁,T^v ₂,T^v ₃,T^v ₄}, we need to find the intersection points k^v along the edges of the tetrahedron. Using an element G_s∈Γ_ℓ, with its center the closest to the center of the element H_v, and using the information from the mesh mapping function we get the spherical coordinates (ϑ ₀,φ ₀) for the node S ₀ of the element G_s. The intersection point k_v of the parametrically defined surface with the edge of the tetrahedron between the nodes T^v ₁ and node T^v ₂ by definition is on that edge, so that:

(T_{2}^{ν} - k_{ν}) \times (T_{2}^{ν} - T_{1}^{ν}) = 0

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 ${(ϑ}_{k_{ν}}, φ_{k_{ν}})$ that correspond to this intersection point according to (9)

k_{ν} = (\begin{array}{c} x_{ν} \\ y_{ν} \\ z_{ν} \end{array}) = (\begin{array}{c} \sum_{l = 0}^{W} \sum_{m = - l}^{l} C_{l, x}^{ℓ, m} {\tilde{Y}}_{l}^{m} (ϑ_{k_{ν}}, {φ_{k}}_{ν}) \\ \sum_{l = 0}^{W} \sum_{m = - l}^{l} C_{l, y}^{ℓ, m} {\tilde{Y}}_{l}^{m} ({ϑ_{k}}_{ν}, {φ_{k}}_{ν}) \\ \sum_{l = 0}^{W} \sum_{m = - l}^{l} C_{l, z}^{ℓ, m} {\tilde{Y}}_{l}^{m} ({ϑ_{k}}_{ν}, {φ_{k}}_{ν}) \end{array}) .

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 (ϑ ₀,φ ₀). 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 ₁,T^v ₂,T^v ₃,T^v ₄} 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 ₁,k^v ₂,k^v ₃}, 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 ₁,k^v ₂,k^v ₃,k^v ₄}, 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).

Fig. 1. Demonstration of the tetrahedral element subdivision scheme used for the two different cases of surface-tetrahedron intersection. (a) In the case when one node denoted (+) is on the one side of the surface while the three others denoted by (-) on the other partition created by the surface. Having found the three intersection points along the adjacent edges of the node, denoted by (x) we divide the tetraherdon into eight subelements shown exploded in (c). (b) In this case, two nodes lay on each side of the surface dividing the tetrahedron into two pentahedral prisms, one on each side of the surface. Having found the four intersection points k^v ₁…k^v ₄, denoted (x) along the edges connecting each the nodes on the one side (+) of the surface to those on the other (-), we cut the element into six tetrahedra subelements shown exploded in (d).

Download Full Size | PDF

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,

Ξ (γ) = {∥ g - A (γ, μ_{a}, D) ∥}_{2}^{2} = {∥ g - A (p) ∥}_{2}^{2}

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

find γ^{*} so that : Ξ (γ^{*}) = min {∥ g - A (γ^{*}, μ_{a}, D) ∥}_{2}^{2}

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

γ^{n + 1} = γ^{n} + {(J_{n}^{T} J_{n} + \land)}^{- 1} J_{n}^{T} (g - A (γ^{n}, μ_{a}, D)) .

where Λ, is a regularisation matrix.

4.1. Jacobian calculation

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

J = \frac{\partial A (γ, μ_{a}, D)}{\partial γ},

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

S (γ, μ_{a}, D) ≔ K (γ, D) + C (γ, μ_{a}) + \frac{1}{2} α E + iω B

we obtain:

S \frac{\partial Φ}{\partial γ} = - \frac{\partial S}{\partial γ} Φ

which leads to

\frac{\partial Φ}{\partial γ} = - S^{- 1} \frac{\partial S}{\partial γ} Φ .

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

J = M [\frac{\partial Φ}{\partial γ}] = - M S^{- 1} \frac{\partial S}{\partial γ} Φ = - Φ^{*} \frac{\partial S}{\partial γ} Φ

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 γ_kare 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 γ_k are derived in equations (22)–(32)of [26].

{(\frac{\partial S}{{\partial γ}_{k}})}_{ij} = (δ μ_{a}) \sum_{p, p^{'}} \int_{Γ} γ_{k} υ_{p} \cdot υ_{p}^{'} + (δ D) \sum_{p, p^{'}} \int_{Γ} γ_{k} \nabla υ_{p} \cdot \nabla υ_{p}^{'}

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:

J_{sd, k} = - ({δμ}_{a}) Φ_{Γ}^{* (d)} G^{k} Φ_{Γ}^{(s)} - (δD) Φ_{Γ}^{* (d)} F^{k} Φ_{Γ}^{(s)}

where

G_{p p^{'}}^{k} = \int_{Γ} γ_{k} φ_{p} (r) \cdot φ_{p'} (r)

F_{pp'}^{k} = \int_{Γ} γ_{k} \nabla φ_{p} (r) \cdot \nabla φ_{p'} (r)

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_p of the volume discretisation where the node S_s of the surface belongs to, and recall the shape functions of this element.

Φ (S_{s}) = \sum_{l = 1}^{4} = Φ (l) υ_{l} (S_{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:

\tilde{g} = (\begin{array}{c} Re (\ln (g)) \\ Im (\ln (g)) \end{array})

which reflects to the scaled Jacobian J̃ defined as

\tilde{J} = (\begin{array}{c} Re [diag (\frac{1}{A (γ, μ_{a}, D)}) J (γ)] \\ Im [diag (\frac{1}{A (γ, μ_{a}, D)}) J (γ)] \end{array})

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₂ 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 (µ_a,2µ′_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

\tilde{Ξ} (γ) = {∥ g^{meas} - A (γ, μ_{a}, D) - g^{ref} + A (p^{ref}) ∥}_{2}^{2}

Where A(p ^ref) represents the solution of the numerical model without any blobs inside and homogeneous optical properties of µ _a=0.0078mm ^-1 and µ′_s=1.065mm ^-1

Fig. 2. The experimental setup used. The sources are denoted by “x” and detectors are de- noted “o”. We also have drawn the positions of the two blobs the one on the left with double the absorption of the background and the one on theright that has double the scattering.

Download Full Size | PDF

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

{‖ g - A (γ^{n}) ‖}_{2}^{2} / | {‖ g - A (γ^{1}) ‖}_{2}^{2}

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.

Fig. 3. Comparison between the target, the reconstruction with the pixel based and the shape based method in vertical cross sections along the x-z plane (a) The real position of the µ _a target in a slice along the y=0 plane and (b) the real position of the µ′s target in a slice along the y=-12.3. (c) The y=0 slice of the recovered µ _a image from the use of a pixel based method as described in [34].(d) The y=-12.3 slice of the recovered µ′_s image from the use of a pixel based method as described in [34]. We notice that the pixel based reconstruction is successful but the contrast is smoothed out due to the regularisation used and the homogeneous background contains artifacts.(e) A view along the y-axis of the recovered µ _a object from the shape based reconstruction in solid color and the initial guess for the object as an outline. (f) A view along the y-axis of the recovered µ′_s object from the shape based reconstruction in solid color and the initial guess for the object as an outline. We notice that the recovered shapes and locations of the object are very close to the expected ones.

Download Full Size | PDF

Fig. 4. Comparison of the view along the central x-y plane between target, the reconstruction with the pixel based and the shape based method As in figure 3 for horizontal slices along the z=0 plane, the real positions of the targets in (a) and (b) and the result from the pixel based method in (c) and (d), and a view along the z-axis (e) and (f) for the shape based method.

Download Full Size | PDF

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 ₁ and Γ^target ₂ 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 Γⁿ_p at iteration n and is defined as

$H aus (Γ_{p}^{t \arg et}, Γ_{p}^{n}) = \max_{v \in Γ_{p}^{t a r g e t}} \max_{v' \in Γ_{n}^{p}} ∣ v - v' ∣$

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.

C (Γ_{p}^{t a r g e r}, Γ_{p}^{n})=| Γ_{p c}^{t a r g e r} - Γ_{p c}^{n} |

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.

Fig. 5. Metrics used to assess the quality of the reconstruction, plotted for the iterations of the shape evolution. (a) The normalized residual ‖g-A(γⁿ)‖² ₂/|g-A(γ ¹)‖² ₂. (b) the volume ratios between targets and reconstructed shapes. (c) the Hausdorff measure for the mismatch between the surfaces of the shapes and (d) the distance between the center of the target objects and the reconstructed ones.

Download Full Size | PDF

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.

Fig. 6. Reconstructions from experimental data, using different initial positions, as presented in the media files attached to the paper. The first row of the graph consists of the initial guesses in yellow and the experimental targets, red for absorption and blue for scattering displayed in Fig. 6(a) (media 1) Fig. 6(b) (media 2) and Fig. 6(c) (media 3). In the second row we have the reconstruction that resulted from those different initial guesses respectively, with the targets recovered successfully in each case.

Download Full Size | PDF

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]

3D shape based reconstruction of experimental data in Diffuse Optical Tomography

Abstract

1. Introduction

2. Forward model

3. Discretization with FEM

3.1. Classification of FEM mesh nodes

3.2. Classification of mesh elements

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

3.4. Subdivision of the tetrahedron in subtetrahedra

3.5. Computation of the system matrix

4. Inverse Problem

4.1. Jacobian calculation

4.2. Scalings

5. Results

6. Conclusion

Acknowledgments

References and links

Supplementary Material (3)

Cited By

Figures (6)

Equations (38)

Optics Express