Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Med Imaging. Author manuscript; available in PMC 2008 December 1.
Published in final edited form as:
PMCID: PMC2591024

Estimation and Statistical Bounds for Three-Dimensional Polar Shapes in Diffuse Optical Tomography

Gregory Boverman, Member, IEEE, Eric L. Miller, Senior Member, IEEE, Dana H. Brooks, Senior Member, IEEE, Qianqian Fang, Member, IEEE, and David A. Boas


Voxel-based reconstructions in Diffuse Optical Tomography (DOT) tend to produce very smooth images due to the attenuation of information of high spatial frequency. This then causes difficulty in estimating the spatial extent and contrast of anomalous regions such as tumors. Given an assumption that the target image is piecewise constant, we can employ a parametric model to estimate the boundaries and contrast of an inhomogeneity directly. In this paper, we describe a method to directly reconstruct such a shape boundary from diffuse optical measurements. We parameterized the object boundary using a spherical harmonic basis, and derived a new method to compute sensitivities of measurements with respect to shape parameters, based on an integral equation formulation of the forward problem. We introduced a centroid constraint to ensure uniqueness of the combined shape/center parameter estimate, and a projected Newton method was utilized to optimize the object center position and shape parameters simultaneously. Using the shape Jacobian, we also computed the Cramér-Rao lower bound on the theoretical estimator accuracy given a particular measurement configuration, object shape, and level of measurement noise. Knowledge of the shape sensitivity matrix and of the measurement noise variance allows us to visualize the shape uncertainty region in three dimensions, giving a confidence region for our shape estimate. We have implemented our shape reconstruction method, using a finite-difference-based forward model to compute the forward and adjoint fields. Reconstruction results are shown for a number of simulated target shapes, and we investigate the problem of model order selection using realistic levels of measurement noise.

I. Introduction

Diffuse Optical Tomography (DOT) is an emerging imaging modality which is beginning to show great promise, primarily due to its ability to monitor the body's hemodynamic properties in vivo using near-infrared light. Among DOT's advantages over conventional imaging modalities in widespread use are its relatively low cost and the employment of safe, non-ionizing radiation. However, image reconstruction for DOT is difficult and relatively computationally expensive, primarily due to the ill-posedness and nonlinearity of the inverse problem. The ill-posedness is caused by the multiple scattering that photons experience in their path from source to detector, with the result that a typical measurement is affected by the absorption over a large volume of the imaging geometry. If it is assumed that the forward problem can be linearized, it is possible to compute the analytical solution to the regularized inverse problem for a number of specialized imaging geometries [1], but it is not currently possibly to do so for more complex geometries such as the breast. In the case of the fully nonlinear inverse problem, which we address in this paper, analytical methods may be applicable for the solution of the inverse problem [2], but again only for specific imaging geometries, and they require the assumption that the scattering series can be truncated to a particular order. Thus, for realistic imaging geometries, nonlinear optimization methods must be used.

In this paper, we will focus on the application of DOT for breast tumor detection and characterization. Recent experimental work has imaged tumors with DOT and has suggested that tumors may differ sharply from normal tissue in their blood volume, oxygenation, and, perhaps in their scattering properties [3] [4] [5].

The most common approach to the DOT inverse problem in the literature is to estimate an absorption distribution which minimizes a cost functional. This cost functional typically includes a term which represents the mismatch between the actual measurements and the measurements which would have resulted from the estimated absorption distribution, generally computed by means of a forward solver. An additional term is included for the purpose of regularization, penalizing some function of the estimate itself, such as its norm, or the norm of its gradient. The estimation problem, having been reduced to a nonlinear regression problem, can be solved by means of gradient or Newton-based approaches, among others. The functional gradient of the cost function with respect to absorption perturbation can be computed efficiently [6]. The resulting solution is highly dependent on the choice of a regularization parameter, the selection of which has been quite exhaustively discussed in the inverse problems literature [7] [8].

However, using standard approaches, recovery of edges in the DOT inverse problem is extremely difficult. The problem is that the propagation of diffuse light highly attenuates absorption distributions with high spatial frequency, and thus inversion methods must make use of regularization to prevent the inversion result from being overwhelmed by noise. Thus, the preponderance of DOT image reconstructions found in the have little discernible shape [9], [3]. One approach taken in the literature on image reconstruction in order to deal with this this problem is to incorporate total-variation minimization [10] or edge-preserving regularization [11] into the image reconstruction. This approach has also been useful in the related applications of image restoration [12], [13] and image segmentation [14], [15].

In this paper, we will show that, given the assumption that the medium is piecewise constant, we can accurately reconstruct the shape, optical absorption contrast, and location of an anomaly, making use of algorithms with reasonable computational complexity. The shape estimation problem is of great importance in itself and because the accuracy in estimating an object's shape is intimately connected with the accuracy in determining the contrast. Quantitatively accurate contrast estimates may be necessary if the image reconstructions produced by DOT are to be clinically informative.

We recover high spatial-frequency edge information which is attenuated, but not completely lost, by explicitly postulating the existence of a boundary and estimating the shape of this boundary using a nonlinear optimization method. While it is not exactly true that the body, or the breast in particular, is a piecewise-constant medium, it is likely to be composed of distinct regions (i.e. adipose tissue, glandular tissue, tumor tissue), with clear boundaries, and where the mean of the optical properties within each region differs significantly from that of the other regions [16]. In this paper, we consider a simpler problem in which we are concerned only with the estimation of a single inclusion embedded in an otherwise homogeneous medium. We demonstrate the properties of this approach on simulated data with varying levels of measurement noise.

Recent shape-based work has made use of level-set functions [17], active contours [18], and various boundary parametrizations. In two dimensions, Fourier Descriptors [19], [20] and B-Spline models [21] have been employed to parameterize the boundary. For DOT, the first work in this area examined the estimation of the location, contrast and orientation of ellipsoids [22]. More recently, a spherical harmonic shape parametrization has been used in conjunction with the boundary element method (BEM) [23]. In our work, we use spherical harmonics to model fairly intricate three-dimensional polar shapes, as described by Li [15], (a polar shape is one that can be modeled as a single-valued function in spherical coordinates), with respect to some center position. We show that we can efficiently solve the inverse problem, using shape derivative theory and the adjoint method to compute the derivatives of our measurements as we vary our shape parameters. We also introduce a centroid constraint to allow simultaneous estimation of object position and shape.

We will make use of the domain, or the Fréchet derivative of the scattered field with respect to infinitesimal perturbations of the object boundary, as described by Sokolowski and Zolesio [24]. The domain derivative approach was applied to the electromagnetic inverse problem for the case of a perfectly conducting obstacle by Kirsch [25]. Hettlich [26] extended Kirsch's variational approach to compute the Fréchet derivative of the scattered field with respect to boundary perturbations for a number of object boundary conditions. The domain derivative can be computed as the solution of a forward problem with particular boundary conditions defined on the object boundary. Potthast [27] gave an alternative derivation for the scattered field domain derivative, assuming a C2 continuously differentiable boundary. Hohage and Schormann [28] extended this method to compute the Fréchet derivative of the scattered field for perturbation of the boundary of a penetrable obstacle. Our derivation will follow that of Potthast, modified for the case of a penetrable obstacle, and making use of volume integrals rather than surface integrals. We will make use of the adjoint approach [29], [30], [31] in order to efficiently compute the Fréchet derivative of boundary measurements without the necessity of a large number of forward solves.

In addition, knowledge of the shape and center Jacobian matrices allows us to compute the constrained Cramér-Rao lower bound on shape estimation accuracy and to visualize this bound. This bound tells us the best reliability that is possible in our shape reconstruction, given a particular measurement geometry and measurement noise level. This knowledge, apart from giving us a shape “confidence interval”, can help us in our design of instruments and source/detector layouts, if we have prior knowledge that a particular shape reconstruction reliability is desired. In contrast to previous work on shape estimation bounds in image reconstruction problems [21], [19], [32], we derive and implement these bounds specifically for the DOT inverse problem, and for the very challenging problem of fully three-dimensional nonlinear shape estimation. We visualize these bounds using a method introduced by Kirsch [25], extended to the case of three-dimensional polar shapes. As we have introduced a centroid constraint to ensure uniqueness of our object shape/center parametric representation, we make use of the constrained Cramér-Rao lower bound [33] to give us performance bounds on simultaneous estimation of object shape and center.

In our experimental work, we show in simulations that fairly complex objects can be estimated with a high degree of accuracy using DOT data and a reasonable level of noise. We also show shape reconstruction performance where there is a trade-off between model complexity and estimation accuracy, since overly complex models cannot be accurately estimated using a limited set of noisy measurements. We examine the question of model order selection, showing that a sufficient level of model complexity is necessary in order to estimate an absorber's size and contrast accurately, but that the model complexity need not be identical to the object's complexity in order to achieve reasonably accurate estimates.

Finally, we look to the question of shape reconstruction lower bounds, visualizing the shape bounds derived from the parametric Cramér-Rao bounds, using both the unconstrained bound and the constrained bound, where a centroid constraint has been introduced.

II. Shape Parameterization

In order to reduce the complexity of the three-dimensional shape estimation problem, we make the assumption that we are reconstructing polar shapes, which can be expressed as single-valued functions in spherical coordinates, with respect to some center location. In this context, spherical harmonics are a very useful choice for shape parametrization, as they define an orthonormal basis on the unit sphere. Since spherical harmonics are complex functions, we make use of their real and imaginary parts, the Tesseral harmonics, in order to ensure that our shape estimate is strictly real. The Tesseral harmonics of order 4 and below are shown in Fig. 1. We note that the shape representation in [23] is more flexible than that used in our work in that it is not limited to star-shaped objects, but this representation introduces considerable difficulties in the visualization (which we plan to address in a future work).

Fig. 1
Spherical Harmonics of up to fourth order. (a) Real part. (b) Imaginary part.

The spherical harmonics are defined as follows [34]:

equation M1

where θ [set membership] [0, π) and [var phi] [set membership] [0, 2π) and equation M2 is the associated Legendre function. Therefore, the radius with respect to a given center position can be represented as follows [15]:

equation M3

where equation M4 and equation M5. We also note that equation M6 and equation M7. The square-root is employed to assure that the radius will be always non-negative, and a small constant, [sm epsilon], is added to ensure the differentiability of the radial function with respect to the shape parameters. It is useful to note that this shape representation is not unique, since the same shape can be represented using a different center position and appropriately modified shape parameters.

In matrix form, we have the following parametrization of the shape boundary:

equation M8

where the rows of B are made up of our shape basis functions and equation M9.

In the results which follow, it will be necessary to compute the the normal vector at each point of our parametrized shape's surface. An analytical expression for this function can be derived by transformation to rectangular coordinates:

equation M10

If we represent each point on the surface of the shape as f = (x, y, z), then the normal vector, n^, can be computed as:

equation M11


equation M12
equation M13

In order to compute equation M14, it is useful to utilize the following identities:

equation M15

A. Centroid constraint

The shape parametrization described above is not unique, as a given shape S can be described with respect to an arbitrary center point within its interior. To ensure uniqueness, we constrain our shape representation such that the center with respect to which it is defined is identical to the shape's centroid, which is unique and well-defined.

Thus, we impose the following constraints, which can be derived by integrating the expressions for the x, y, and z coordinates of the estimated shape's centroid in spherical coordinates:

equation M16
equation M17
equation M18

III. Forward Modeling

In a given DOT experiment, the medium is illuminated using modulated light by a point source, located at rs and measurements are made at discrete detector positions, with detectors approximated as delta functions at the detector positions. The experiment is then repeated using a point source located at a different position. In what follows, we assume that all measurements are made in the frequency domain, using a single modulation frequency. We assume that photon propagation is well approximated by the diffusion equation [35] and make use of a zero partial-flux Robin-type boundary condition [36]. We also assume that the medium is piecewise constant and consists of two regions, Ωi and Ωe, where Ωi is simply connected, separated by an internal boundary Γi. We denote the exterior boundary by Γe. The photon intensity, u(r) can be modeled as the the solution of the following system of coupled partial differential equations:

equation M19

with the following boundary conditions on Γi:

equation M20

where Di and De are the diffusion coefficients in Ωi and Ωe respectively, αi,e are the absorption coefficients, the modulation frequency is ω, υ is the speed of light in tissue, and the unit outward normal vectors are n^i,e(r). The coefficient β is dependent on the indices of refraction at the exterior boundary. In the remainder of this paper, we will assume that Di = De.

In addition to the above formulation in terms of partial differential equations, us (r), the photon fluence resulting from illumination at a given source position can be formulated as a Fredholm integral equation of the second kind, as follows:

equation M21

where G(r, r′) is the Green's function in the absence of the absorber and Ω = Ωi [union or logical sum] Ωe. In operator notation, we can then express this equation as:

equation M22


equation M23
equation M24

and α = αiαe. For future use, we define the surface integral operator, GΓi as:

equation M25

IV. Optimization and Sensitivity Calculations

In a typical DOT experiment, the medium is illuminated sequentially at each of M source positions, with measurements made at N detector positions. If we arrange our measurements as a vector, the goal of our inverse problem is to minimize the following cost function:

equation M26
equation M27
equation M28

where |yi,j| and Arg(yi,j) are the magnitude and phase, respectively, of the measurement of the fluence due to source i, measured at detector location j. |hi,j| and Arg(hi,j) are the corresponding hypothesized measurements of magnitude and phase, given a forward solver and an estimate of the object shape parameters, p. The covariance matrix of the noise, assumed to be Gaussian, is given by Σn.

In this paper, p has the following form:

equation M29

where cx, cy, cz are the x, y, and z components of the object center position, respectively, α is the absorption contrast, and a is the shape parameter vector.

In the results that follow, we show how to explicitly compute the derivative of the cost function with respect to an object's shape parameters and its center position. Although we make use of spherical harmonic basis functions, the following analysis applies to any representation of polar three-dimensional shapes, and, in fact, the analysis of the immediately following section applies to an arbitrary shape model which is not necessarily polar.

A. The Shape Derivative

In this section, our goal is to derive the shape derivative of the object boundary, which is the change in the total field us(r) with respect to infinitesimal variations in Ωi. We will consider two cases. In the first case, we will give an expression for the shape derivative where Ωi is perturbed by an arbitrary, smooth vector field, V. In practice, we will make use of this approach in computing the sensitivity of our shape represenation to object center position. In the second case, we will restrict ourselves to the case where Γi is polar and parametric. In order to accomplish this objective, we need a few basic ideas from nonlinear functional analysis [37]:

An operator A: UX mapping an open subset of a normed space to a Banach space is Fréchet differentiable if:

equation M30

where A1(h) = o(| |h| |) and equation M31 is a bounded linear operator, the Fréchet derivative of A(r0).

Assuming that A is invertible and continuous in some neighborhood of y0, then it can be shown that [27]:

equation M32

If V is a C2 differentiable vector field defined in some neighborhood of Γi, then we can define the shape derivative of the field due to a particular point source with respect to an infinitesimal perturbation of Γi in direction V as follows:

equation M33

It can be shown that [24], if y = GΩix then the shape derivative of y with respect to a perturbation of Ωi by V is:

equation M34

Using Eqs. 24, 26, and 15, we can compute the shape derivative of the photon fluence as:

equation M35
equation M36

This result shows that computation of the Fréchet derivative is equivalent to solving the forward problem with particular boundary conditions on Γi [26]. Each evaluation of (I + αGΩi)−1y for a given function y requires a forward model computation.

Making use of this expression, we can compute the sensitivity of the measurements to perturbations in the object's center position. Perturbing the x coordinate of our object's center is equivalent to shape transformation by a vector field where V = [ 1 0 0 ] everywhere. Thus the derivative of a measurement with respect to a change in the x component of the center position can be computed as:

equation M37

where n^x(s) is the x component of n^i(s) at s. The sensitivities to changes in the y and z components of the absorber's center position can be computed in an analogous fashion.

We now consider the case where the boundary can be parametrized using polar coordinates and that it is dependent on a vector of discrete parameters, Γi = Γi[r0, θ, [var phi], R(θ, [var phi], a)], where θ [set membership] [0, π) and [var phi] [set membership] [0, 2π). In this case, we can compute the Fréchet derivative of u with respect to changes in the parameter vector, a:

equation M38

The Fréchet derivative of us with respect to the contrast, α, can similarly be computed:

equation M39

B. Adjoint formulation

Unfortunately, actually computing the sensitivity of us with respect to a using Eq. 30 would require one forward solve for each element of a and an additional forward solve for the contrast sensitivity thereby providing the sensitivity information over all points in space. For the purposes of parameter estimation, however, we do not require all of this information. Rather, the optimization approach which we will use requires the derivative of the photon fluence only at the detector locations. Using the adjoint approach [29], [31], then, we can significantly reduce the computational demands of our shape sensitivity calculations.

In what follows, we make use of the following inner product defined over Ω:

equation M40

For the diffusion operator, it is well known that G(r, r′) = G*(r′, r) which leads to the observation that < Gu, υ >Ω = < u,Gυ >Ω. Formally, the adjoint method proceeds by noting that the measurement at detector d due to source s can be written as:

equation M41

Thus, where in the previous section we computed the Fréchet derivative of the field usi excited by source si and measured at all points in space, here we are instead computing the derivative of a scalar measurement, excited at a particular source position, and measured at a single point in space. The measurement operator applied to the forward model solution is:

equation M42

Using Eq. 28, we can compute the shape derivative of this expression as follows:

equation M43

Finally, we show that the shape semiderivative can be computed as the following relatively simple surface integral:

equation M44

where ũd(s) is the solution of the forward problem using a source at the detector position. We show this by writing Eq. 35 in integral form, where equation M45:

equation M46
equation M47

Eq. 38 results from the observation that K(r′, r) = K*(r, r′). Eq. 36 results by defining ũd(r) as:

equation M48

We observe that ũd(r) can be computed using a forward solve, with an “adjoint” source at the detector position.

In the same way, we can manipulate Equations 30 and 31 to compute the sensitivity of discrete boundary measurements of the fluence with respect to perturbations in the shape parameters of a polar, parameterized shape:

equation M49
equation M50

We see that once the forward solutions have been computed for all sources and the adjoint solutions have been computed for all detectors, the above sensitivity calculations can be computed with relatively little computational cost using surface or volume integrals over Ωi or Γi, respectively.

Making use of Eq. 29, extending the result to perturbations in all components of the object's center position, and making use of the adjoint procedure described above, we can derive the following expression for the gradient of a measurement with respect to center position:

equation M51

Again, this formula is an integral over the surface of the object, and is valid for computing the derivative of our measurements with respect to changes in an absorbing object's center position, for arbitrarily complicated shapes.

C. Optimization

In the shape optimization, we make use of a projected Newton method [38], which is the projection of the Levenberg-Marquardt method onto the manifold defined by Eqs. 9, 10, and 11. Specifically, we employ the following algorithm:

Algorithm 1

  1. Choose values for c0, the initial center position, α0, the initial contrast, and a0, the initial shape parameter vector.
  2. Compute the source and detector solutions for the current shape estimate: ci, αi, ai
  3. Compute the cost function: equation M52
  4. Compute Ji = [ Jc,iJα,iJa,i ]
  5. Update the estimate:
    equation M53
    where we choose λ and T to minimize the cost function of the new estimate.
  6. Return to step 2 until convergence.

In Eq. 43, U(t) is an orthonormal basis for the null space of the gradient matrix of the constraint set defined by Eqs. 9, 10, 11, and t is the length of the step taken in the search direction. The Jacobian matrices, Jc,i, Jα,i, and Ja,i are obtained by computing Eqs. 42, 41, and 40, respectively, for all source-detector measurement pairs. We note that we assume that y and h are comprised of measurements of amplitude and phase and we make use of the chain rule in order to compute the Jacobian matrices with respect to these measurement types. We approximate the integral in Eq. 43 numerically. This optimization is an extension of the standard Levenberg-Marquardt method, where the search direction is projected onto a lower-dimensional manifold defined by the constraint set.

For the purpose of comparison, we also make use of an algorithm in which the object center position is known a priori. In this case, we use an unconstrained Levenberg-Marquardt algorithm and optimize only the shape parameters and the object contrast:

Algorithm 2

  1. Choose values for α0, the initial contrast, and a0, the initial shape parameter vector.
  2. Compute the source and detector solutions for the current shape estimate: αi, ai
  3. Compute the cost function: equation M54
  4. Compute Ji = [ Jα,iJa,i ]
  5. Update the estimate:
    equation M55
    where we choose λ and T to minimize the cost function of the new estimate.
  6. Return to step 2 until convergence.

V. Voxel-Based Nonlinear Reconstruction

We have also implemented a more standard image reconstruction algorithm in which we estimate the absorption at each point in space by assuming that it can be decomposed as a linear combination of piecewise-constant basis functions:

equation M56

where χn (r) is the characteristic function of voxel n.

Estimating the absorption at each point in space, given this assumption, is then equivalent to estimating the vector v = [ υ1υ2υn ]T. It can be shown that the Fréchet derivative of a measurement hs,d with respect to a small perturbation in v is [6]:

equation M57

where, again, us(r) and ũd(r) are the forward and adjoint solutions for source s and detector d, respectively.

We then make use of a nonlinear optimization algorithm to estimate v, with Tikhonov regularization applied to stabilize the inversion. Our optimization algortithm, then, strives to minimize the following quadratic cost function:

equation M58

where λ is the regularization parameter and L is the finite-difference approximation of the Laplacian, used to introduce smoothing of p. We implement L as a sparse matrix. Specifically, we follow the procedure listed below:

Algorithm 3

  1. Choose a value for v0, the initial guess. As our initial guess, we assume that v = υ[ 1 1 … 1 ]T, where υ is known a priori.
  2. Compute the source and detector solutions for the current value of v.
  3. Compute the cost function: equation M59
  4. Compute the Jacobian matrix Jv,i, using 46 for all source-detector measurement pairs.
  5. Update the estimate:
    equation M60
    where we choose T to minimize the cost function of the new estimate.
  6. Return to step 2 until convergence.

VI. Cramér-Rao Lower Bound and Shape Bounds

An advantage of being able to explicitly compute shape sensitivities is that we can compute a lower bound on shape parameter estimation accuracy. For an unbiased estimator, we can compute a lower bound on the estimation covariance using the Fisher Information Matrix (FIM):

equation M61

where the FIM is:

equation M62

and f(y, p) is the joint probability density of the observation vector y and the parameter vector p.

If our estimator is asymptotically Gaussian, we have the following simplified form for the FIM for the estimation of a:

equation M63

where Jp is the Jacobian of our measurements with respect to p, each element of which can be computed using Eqs. 40, 41, and 42.

In our case, where there are additional constraints on p, it is necessary to use the constrained Cramér-Rao bound. If c(p) is a set of differentiable constraints and Jc(p) is its gradient matrix, then it can be shown that [33]:

equation M64

where the constrained Fisher information matrix is:

equation M65

and U is a matrix whose columns form an orthonormal basis for the null-space of Jc(p).

Given a lower bound for the parameter covariance matrix, Ep[(p^p)(p^p)T] we wish to determine bounds on the shape itself. To accomplish this end, for each point on the boundary, we compute the maximum deviation such that the parameter vector is within a given confidence interval. Thus, we solve the following constrained optimization problem for each point on the boundary:

equation M66

where γ is our confidence interval. In this equation, the 0 block, with dimension 1 × 4 is used to select only those shape parameters used to encode the object radius, ignoring those which relate to object position and contrast. This expression can be used to computed the maximal perturbation in the object radius, where [p with tilde] = p^p. The same procedure can be used to give us a confidence interval on the object center position and on the contrast.

VII. Results and Discussion

In Fig. 2 (a), we show a shape generated using spherical harmonic coefficients of order 3 and below, for which we have simulated diffuse optical measurements. The object is embedded in a cubical region filled with a diffusing medium 6 cm on a side, and we make use of a transmission geometry, with 25 sources and 25 detectors uniformly distributed on the z = 0 cm and z = 6 cm planes, respectively. A background absorption of 0.05 cm−1 was assumed, and we assumed a (fairly high) absorption contrast of 0.15 cm−1 in order that nonlinear effects would be significant in the image reconstruction problem. The true absorber is centered at (3 cm, 3 cm, 3 cm). We initialize our reconstructions with a spherical object 1.6 cm in diameter centered at (4 cm, 4 cm, 4 cm), with an absorption contrast of 0.07 cm−1.

Fig. 2
Randomly generated third-order shape and performance of a voxel-wise algorithm in reconstructing this shape. (a) True third-order absorbing object. (b) 50% absorption isosurface for the voxel-wise reconstruction. (c) Convergence of the cost function for ...

In all of the examples described in this paper, we make use of a finite-difference-based forward model, discretized using a uniform grid, with voxels 0.2 cm on a side. Partial-volume effects are taken into account by computing the volume of the intersection of the object with each voxel, making use of integration in spherical coordinates. To compute the forward solution, we discretize the PDE described by Eqs. 12 as a large, sparse matrix, and solve using the GMRES iterative method [39].

For the purpose of comparison, we show in Figs. 2 (b) and 2 (c) a voxel-based nonlinear reconstruction of the absorption, computed using Algorithm 3, for data produced by the absorbing object shown in Fig. 2 (a), with amplitude-dependent Gaussian noise added such that the SNR of the amplitude measurements was 80 dB, and Gaussian noise with a standard deviation of 0.1° was added to the phase measurements. The regularization parameter was estimated by means of the L-curve method [40], [7], [41]. We note that the peak absorption contrast is estimated fairly accurately, but, as is generally the case with DOT reconstructions, the reconstructed absorber is quite smooth and is lacking distinct edges. Thus, if the reconstructed image in Fig. 2 (b) represented a tumor, it would be very difficult to definitively estimate the tumor's volume, or for that matter, its actual absorption contrast.

In Fig. 2 (d) we show the convergence behavior of our voxel-wise reconstruction procedure, as a function of iteration, noting that the cost function which we will reference in the remainder of this paper is the square-root of Eq. 19. Fig. 2 (b) shows the 50% isosurface of the reconstructed absorber. We notice that the position, orientation, and general shape of the object are accurately reconstructed, but the fine details of the shape are missed. Reducing the regularization parameter may help to reconstruct these fine details, but the level of image artifacts worsens.

In contrast, the results of the shape-based image reconstruction, using the same data as for the voxel-wise reconstruction, are shown in Fig. 3, which illustrates the evolution of our shape estimation in the first five iterations of the algorithm. We note that, for this case, at least visually, the true shape, position, and orientation of the absorber are reconstructed quite accurately after five iterations of the algorithm. As illustrated in Fig. 3 (d), the cost function continued to decrease further until approximately the twenty-fifth iteration. Fig. 3 (e) shows that the contrast of the absorber is estimated to a very high degree of accuracy after five iterations, but that the convergence rate for the contrast seems to drop somewhat at this point, requiring approximately twenty more iterations for the contrast to converge to its true value. This is typical in our simulations; there seems to be a subtle interplay between accurate estimation of shape and estimation of contrast.

Fig. 3
Evolution of third-order shape estimate and performance of shape-based reconstruction algorithm: (a) Iteration 1. (b) Iteration 3. (c) Iteration 5. (d) Convergence of the cost function. (e) Convergence of the absorption contrast estimate.

Fig. 3 (d) also compares the convergence of the algorithm in the cases where the object center position is known a priori to the case when it is not assumed to be known. In the former case, Algorithm 2 is used, and, in the latter, a projected-Newton method, Algorithm 1, is employed. In both cases, the algorithm converges to a cost considerably lower than the minimum achieved by the voxel-based reconstruction approach. In addition, as expected, the rate of convergence is somewhat slower in the case where we have less a priori information, but the minimum is still reached in approximately 25 iterations. In the case where the center position was known, our estimator was able to reach a value within one standard deviation of the expected Chi-squared error due to noise, meaning that essentially all of the information was extracted from the measurements. This minimum was not reached in the case where the center was not assumed known, and the reason for this is that our projected-Newton algorithm did not satisfy the centroid constraint exactly, and thus the estimated center tended to drift slightly from the true center position as the algorithm progressed (with the shape parameters compensating accordingly).

An important question is whether shape-based image reconstruction methods are robust even in the very realistic circumstance where the medium is not exactly piecewise constant. To explore this question, we have conducted simulations with variation in the medium's background absorption coefficient, where the true absorbing object is shown in Fig. 2 (a). We generated Gaussian independent identically distributed random noise fields with two different values chosen for the noise variance. These random noise fields were added to the absorption throughout the image, including within the absorbing anomaly, with the resulting absorption images for the z = 3.0 cm slice of the medium shown in Figs. 4 (a) and 4 (b). In Figs. 4 (a) and 4 (b) we have simulated random processes with relatively high and low variance, with standard deviations of 0.03 cm−1 and 0.005 cm−1, respectively, where the background absorption was 0.05 cm−1 respectively. The results of shape-based reconstructions, using Algorithm 1, for the levels of background variation in Figs. 4 (a) and 4 (b) are shown in Figs. 4 (c) and 4 (d), respectively. The volume of the actual object, shown in Fig. 2 (a), is 2.94 cm3. For the case of low-variance background variation, the estimated volume was 2.6 cm3 and the estimated absorption contrast was 0.17 cm−1 (the true contrast was 0.15 cm−1), giving errors −11.6% and 13.3% in volume and contrast, respectively. In the case of high-variance background absorption variation, the errors in estimated object volume and contrast were −26.5% and 33.3%, respectively.

Fig. 4
Background variation added to investigate the effect of heterogeneous absorption: z = 3.0 cm slice shown. (a) High variance background variation. (b) Low variance background variation. (c) Reconstruction in case of high variance background variation. ...

For the purpose of comparison, we have computed voxel-based reconstructions of absorption, using Algorithm 3, for the cases shown in Figs. 4 (a) and 4 (b). The isosurfaces representing a value of absorption 50% below the peak perturbation are shown in Figs. 4 (e) and 4 (f), respectively. We note that the voxel-based reconstruction algorithm estimated approximately the correct location, contrast(in both cases, a peak contrast of 0.15 cm−1 was estimated), and orientation for the absorbing object, but, again, the reconstructions were excessively smooth and the true shape of the boundary was not recovered.

Next we study model-order selection while reporting the performance of Algorithm 1 for a more complex shape. We also increased the amplitude and phase noise such that the amplitude SNR was only 40 dB and the standard deviation of the phase noise was 1.0°. Fig. 5 (a) shows a randomly-generated shape, using spherical harmonic coefficients of order 8 and below. This example is intended to approximate the spiculated nature of tumors in the breast. Fig. 5 (b)-(j) shows the estimated shape as we vary the order of the model used in the reconstruction. As in the previous example, the absorption contrast was 0.15 cm−1. We note that we are able to reconstruct the fine details of the absorber shape, although it is not visually clear whether there is reconstruction accuracy improvement for models of order 5 and higher.

Fig. 5
Shape reconstruction as a function of model order (a) True eighth-order absorbing object (b) Order 0 estimate (c) Order 1 (d) Order 2 (e) Order 3 (f) Order 4 (g) Order 5 (h) Order 6 (i) Order 7 (j) Order 8

In Fig. 6, we report a quantitative analysis of the results shown in Fig. 5. In Fig. 6 (a) we see that the cost function continues to decrease as the model-order used in the reconstruction is increased until order 6. For orders 6 and above, the cost is approximately constant. Fig. 6 (b) shows that the image reconstruction error also decreases as the reconstruction model-order increases until order 6, but that the error increases significantly using models of orders 7 and 8. Thus, in high-noise situations, using an excessively complex model can be counter-productive. Figs. 6 (c) and (d) show the estimation of object volume and contrast as a function of model order, with the true volume and contrast, respectively, depicted as dashed lines. By numerical integration, we computed a true inhomogeneity volume of 4.82 cm3 and, using a model of order 6, which minimized the reconstruction error, the estimated volume was 4.43 cm3, giving an error of −8.1%. Likewise, for the order 6 model, the estimated contrast was 0.164 cm−1, an overestimation of the true contrast by 9.3%. Even in the case of a relatively high level of measurement noise, then, the volume and contrast of a high-contrast absorber with a very complex shape can be estimated quite accurately, with errors of less than 10% in object volume and absorption contrast. Fig. 6 (e) shows the interesting result that, regardless of the model complexity used in the reconstruction, the product of the estimated contrast and the estimated object volume is approximately constant, as suggested by the analysis in [42].

Fig. 6
Quantitative investigation of model order (a) Residual at estimate to which shape reconstruction converged, as a function of model order (b) Reconstruction error as a function of model order (c) Estimated shape volume as a function of model order (d) ...

Finally, we visualize Cramér-Rao-based shape recontruction bounds for the object shown in Fig. 2. In Fig. 7, we visualize the shape reconstruction uncertainty region where, for the purposes of visualization, the simulated phase noise has been increased to having a variance of 1°. The bounds are conputed using Eq. 54, with γ = 0.95. In (a), we show the reconstructed shape bounds where the center position and object contrast are known a priori, corresponding to Algorithm 2. In (b) the same results are shown where the center position and object contrast are not known a priori, corresponding to Algorithm 1, and the centroid constraints, Eqs. 9, 10, and 11 are imposed. In this case, the constrained Cramér-Rao bound must be used. We observe that the unknown center position and unknown contrast do increase the maximal shape uncertainty, but to a relatively modest degree, less than 15%.

Fig. 7
Visualization of the shape uncertainty region for a third-order shape, with a phase noise variance of 1° (a) Uncertainty in the case where the center and absorption contrast are known a priori (b) Uncertainty in the case where the center and absorption ...

We note that the confidence regions shown in Fig. 7 may be somewhat pessimistic, in that we are computing the maximum perturbation at each point of the boundary given that the shape parameters are within a 95% confidence interval. It may be that a method, such as that described in [19], and extended to three dimensions, may give a somewhat less conservative estimate of the shape confidence region.

VIII. Conclusion and Future Work

In this paper, we have shown that highly accurate estimation of object boundaries is possible in Diffuse Optical Tomography, an imaging modality with typically very low spatial resolution, even with realistic levels of noise added to the measurements. We have demonstrated a Newton-type image-reconstruction algorithm for the nonlinear reconstruction of absorbing inhomogeneities which are simply connected and polar in shape. With the shape boundaries represented using the real and imaginary parts of the spherical harmonic functions, our algorithm simulataneously estimates object shape, center position, and contrast, making use of constrained estimation with a centroid constraint to ensure the uniqueness of the shape/center representation. We are able to compute the sensitivity of the cost function to perturbations in object shape, center position, and contrast directly, by making use of shape derivative theory, thus obviating the need for computationally expensive finite-difference computations. Given knowledge of these sensitivities, and making use of the constrained and unconstrained Cramér-Rao lower bounds, we are able to compute a lower bound on shape, center, and contrast uncertainty for a given absorbing object, source-detector configuration, and level of measurement noise. In addition to giving us a confidence region on the shape estimates which we obtain, this information can help us to design sensors and to determine the best-case performance of an imaging system with known characteristics.

We hope that this work will be of use in improving the resolution of DOT and in making the results more quantitatively useful, particularly in the estimation of the contrast and size of tumors in the breast. Future work will extend the result in this paper to the problem of estimation of objects using measurements at multiple wavelengths, and incorporating spectroscopic constraints. In addition, we will consider in more detail the problem of simultaneous determination of the shape and contrast of an object in addition to estimation of unknown calibration parameters.


GB, DHB, and ELM acknowledge the support of CenSSIS, the Center for Subsurface Sensing and Imaging System, under the Engineering Research Centers Program of the National Science Foundation (award number EEC-9986821). DAB was supported by the National Institutes of Health. GB thanks Prof. David Isaacson for helpful discussions and advice.


1. Markel VA, Schotland JC. Inverse problem in optical diffusion tomography. I. Fourier-Laplace inversion formulas. J Opt Soc Am A. 2001;18(6):1336–1347. [PubMed]
2. Markel VA, O'Sullivan JA, Schotland JC. Inverse porblem in optical diffusion tomography. IV. Nonlinear inversion formulas. J Opt Soc Am A. 2003;20(5):903–912. [PubMed]
3. Pogue BW, Poplack SP, McBride TO, Wells WA, Osterman KS, Osterberg UL, Paulsen KD. Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: pilot results in the breast. Radiology. 2001;218(1):261–266. [PubMed]
4. Cerussi AE, Berger AJ, Bevilacqua F, Shah N, Jakubowski D, Butler J, Holcombe RF, Tromberg BJ. Sources of absorption and scattering contrast for near-infrared optical mammography. Academic Radiology. 2001;8(3):211–218. [PubMed]
5. McBride TO, Pogue BW, Gerety ED, Poplack SR, Osterberg UL, Paulsen KD. Spectroscopic diffuse optical tomography for the quantitative assessment of hemoglobin concentration and oxygen saturation in breast tissue. Appl Opt. 1999;38(25):5480–5490. [PubMed]
6. Arridge SR. Photon measurement density functions. Part 1: Analytical forms. Applied Optics. 1995;34(31) [PubMed]
7. Hansen PC. Rank-Deficient and Discrete Ill-Posed Problems. Philadelphia, PA: SIAM Press; 1998.
8. Vogel CR. Computational Methods for Inverse Problem. SIAM Press; 2002.
9. Jiang H, Paulsen KD, Osterberg UL, Pogue BW, Patterson MS. Optical image reconstruction using frequency-domain data: simulations and experiments. J Opt Soc Am A. 1996;13(2):253–266.
10. Paulsen KD, Jiang H. Enhanced frequency-domain optical image reconstrution in tissues through total-variation minimization. Applied Optics. 1996;35(19):3447–3458. [PubMed]
11. Charbonnier P, Blanc-Fèraud L, Aubert G, Barlaud M. Deterministic edge-preserving regularization in computed imaging. IEEE Transactions on Image Processing. 1997;6(2):298–311. [PubMed]
12. Vogel CR, Oman ME. Fast, robust total variation-based reconstruction of noisy, blurred images. IEEE Transactions on Image Processing. 1998;7(6):813–824. [PubMed]
13. Li Y, Santosa F. A computational algorithm for minimizing total variation in image restoration. IEEE Transactions on Image Processing. 1996;5(6):987–995. [PubMed]
14. Freedman D, Radke RJ, Zhang T, Lovelock DM, Chen GTY. Model-based segmentation of medical imagery by matching distributions. IEEE Transactions on Medical Imaging. 2005;24(3):281–292. [PubMed]
15. Li J. PhD dissertation. University of Michigan; Ann Arbor, MI: 2002. Three dimensional shape modeling: Segmentation, reconstruction and reconstruction.
16. Brooksby B, Pogue BW, Jiang S, Dehghani H, Srinivasan S, Kogel C, Tosteson TD, Weaver J, Poplack SP, Paulsen K. Imaging breast adipose and fibroglandular tissue molecular signatures by using hybrid MRI-guided near-infrared spectral tomography. Proceedings of the National Academy of Science, USA. 2006;103(23):8828–8833. [PubMed]
17. Santosa F. A level-set approach for inverse problems involving obstacles. ESAIM: Control, Optimization and Calculus of Variations. 1996;1:17–33.
18. Feng H, Karl WC, Castanon DA. A curve-evolution approach to object-based tomographic reconstruction. IEEE Transactions on Image Processing. 2003;12(1) [PubMed]
19. Ye JC, Bresler Y, Moulin P. Cramer-Rao bounds for 2-D target shape estimation in nonlinear inverse scattering problems with application to passive radar. IEEE Trans Ant Prop. 2001;49(5):771–783.
20. Kolehmainen V, Vaukohnen M, Kaipio JP, Arridge SR. Recovery of piecewise constant coefficients in optical diffusion tomography. Optics Express. 2000;7(13):468–480. [PubMed]
21. Hero AO, Piramuthu R, Fessler JA, Titus SR. Minimax emission computed tomography using high-resolution anatomical information and B-spline models. IEEE Transactions on Information Theory. 1999;45(3):920–938.
22. Kilmer ME, Miller EL, Boas D, Brooks D. A shape-based reconstruction technique for DPDW data. Optics Express. 2000;7(13):481–491. [PubMed]
23. Zacharopoulos AD. PhD dissertation. University of London; London, UK: 2005. Three-dimensional shape-based reconstructions in medical imaging.
24. Sokolowski J, Zolesio JP. Introduction to Shape Optimization: Shape Sensitivity Analysis. Berlin: Springer-Verlag; 1992.
25. Kirsch A. The domain derivative and two applications in inverse scattering theory. Inverse Problems. 1993;9(1):81–96.
26. Hettlich F. Frechet derivatives in inverse obstacle scattering. Inverse Problems. 1995;11:371–382.
27. Potthast R. Fréchet differentiability of boundary integral operators in inverse acoustic scattering. Inverse Problems. 1994;10:431–447.
28. Hohage T, Schormann C. A Newton-type method for a transmission problem in inverse scattering. Inverse Problems. 1998;14:1207–1227.
29. Norton SJ. Iterative inverse scattering algorithms: Methods of computing frechet derivatives. J Acoust Soc Am. 1999 November;106(5):2653–2660.
30. Litman A, Lesselier D, Santosa F. Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set. Inverse Problems. 1998;14:685–706.
31. Feijóo GR, Oberai AA, Pinsky PM. An application of shape optimization in the solution of inverse acoustic scattering problems. Inverse Problems. 2004;20:199–228.
32. Ye JC, Bresler Y, Moulin P. Asymptotic global confidence regions in parametric shape estimation problems. IEEE Trans Inf Theory. 2000;46(5):1881–1895.
33. Stoica P, Ng BC. On the Cramér-Rao bound under parametric constraints. IEEE Signal Processing Letters. 1998;5(7):1285–1301.
34. Dennery P, Krzywicki A. Mathematics for Physicists. Mineola, NY: Dover; 1995.
35. Arridge SR. Optical tomography in medical imaging. Inverse Problems. 1999;15:R41–R93.
36. Haskell RC, Svaasand LO, Tsay TT, Feng T, McAdams MS, Tromberg BJ. Boundary conditions for the diffusion equation in radiative transfer. J Opt Soc Am A. 1994 October;11(10):2727–2741. [PubMed]
37. Berger MS. Nonlinearity and Functional Analysis. New York, NY: Academic Press; 1977.
38. Luenberger DG. Optimization by Vector Space Methods. New York, NY: John Wiley & Sons; 1969.
39. Saad Y, Schultz MH. GMRES: A generalized minimal-residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput. 1986;7:856–869.
40. Hansen PC. Analysis of discrete ill-posed problems by means of the l-curve. SIAM Review. 1992;34(4):561–580.
41. Hansen PC, O'Leary DP. Use of the L-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific Computing. 1993;14(6):1487–1503.
42. Boas DA, O'Leary MA, Chance B. Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytical solution and applications. Proceedings of the National Academy of Science, USA. 1994;91:4887–4891. [PubMed]