PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Opt Soc Am A Opt Image Sci Vis. Author manuscript; available in PMC Jun 3, 2012.
Published in final edited form as:
J Opt Soc Am A Opt Image Sci Vis. Mar 1, 2012; 29(3): 321–330.
PMCID: PMC3366196
NIHMSID: NIHMS368989
Information loss and reconstruction in diffuse fluorescence tomography
Petra Bonfert-Taylor,1,2* Frederic Leblond,2 Robert W. Holt,2 Kenneth Tichauer,2 Brian W. Pogue,2 and Edward C. Taylor1,2
1Department of Mathematics, Wesleyan University, 265 Church Street, Middletown, Connecticut 06459, USA
2Thayer School of Engineering, Dartmouth College, 8000 Cummings Hall, Hanover, New Hampshire 03755, USA
*Corresponding author: Petra.B.Taylor/at/dartmouth.edu
This paper is a theoretical exploration of spatial resolution in diffuse fluorescence tomography. It is demonstrated that, given a fixed imaging geometry, one cannot—relative to standard techniques such as Tikhonov regularization and truncated singular value decomposition—improve the spatial resolution of the optical reconstructions via increasing the node density of the mesh considered for modeling light transport. Using techniques from linear algebra, it is shown that, as one increases the number of nodes beyond the number of measurements, information is lost by the forward model. It is demonstrated that this information cannot be recovered using various common reconstruction techniques. Evidence is provided showing that this phenomenon is related to the smoothing properties of the elliptic forward model that is used in the diffusion approximation to light transport in tissue. This argues for reconstruction techniques that are sensitive to boundaries, such as L1-reconstruction and the use of priors, as well as the natural approach of building a measurement geometry that reflects the desired image resolution.
Optical tomography with diffused light is an imaging modality that uses light in the red and near-infrared (NIR) wavelength band to transilluminate biological tissue for the purposes of medical diagnostics (see, e.g., [13]). Signal detection through up to several centimeters of tissue is feasible in that range because scattering dominates over absorption (e.g., from hemoglobin), with the latter being orders of magnitude smaller when compared to other wavelengths in the visible part of the electromagnetic spectrum. Optical tomography is implemented via the mathematical reconstruction of absorption and scattering tissue properties from surface measurements following excitation with light sources such as lasers [4]. Since red and NIR radiation are strongly scattered by biological tissue, one can assume photons are isotropically scattered with-in a few millimeters from the boundary, and as such light transport modeling can be achieved with reasonable accuracy in the diffusion approximation regime of the radiative transport equation [5,6].
Of interest here are two types of imaging modalities, namely, diffuse optical tomography (DOT) and diffuse fluorescence tomography (FT), with the latter often being referred to as fluorescence molecular tomography (FMT). The optical contrast targeted for reconstruction in DOT is typically associated with tissue chromophores such as oxygenated and deoxygenated hemoglobin and sometimes with water and lipids [13]. However, FT deals with the optical contrast that is induced by fluorescent reporters following the systemic injection of a compound targeting specific molecular events associated with pathological processes [7] (see also [8] for a review). It is a well-known fact that the spatial resolution and the contrast recovery associated with optical tomography methods present engineering and mathematical challenges. Fundamental spatial resolution limits in diffused imaging are set by the mean photon free path or the inverse of the transport coefficient, which is typically smaller than 1 mm. However, although several studies have considered the actual spatial resolution and contrast recovery limits observed in re-constructed images, there is no universal criterion that has been derived to quantify them partly because they depend on the particular details of the imaging geometry and the biological application of interest.
A critical aspect impacting the imaging capabilities of optical tomography is the fact that the associated inverse problem is ill posed mathematically. More precisely, diffuse tomography models are highly ill conditioned, implying that even minute discrepancies between model predictions and experimental data might be a dominating factor in non-regularized reconstructions. Typical sources of data–model mismatch include discretization errors, stochastic noise, light transport modeling inaccuracies, and systematic errors. However, while the use of regularization allows one to maximize the amount of physical information recovered from the data, ill-posedness implies that a significant portion of the information that is available in principle can be lost in the process. As demonstrated in, e.g., [911], the amount of information that is lost through regularization can be minimized by improving the light detection geometry by maximizing the level of signal to noise and by using optical data types leading to less ill-posed problems. In these studies, the tool of choice for evaluating the degree of ill-posedness consists of a singular value analysis of the matrices derived from light transport solutions. Typically, the most favorable imaging problems will be those for which the singular value curves (as a function of mode order) have a slope that is as small as possible, implying that high-frequency components (e.g., associated with shot noise) will have less of a tendency to corrupt the images. A conclusion of those studies is that system design optimization should be done ensuring that image reconstruction can be achieved using as many singular modes as possible, especially when the number of available measurements is limited.
The limitations, however, associated with optical tomography methods are also impacted by the fact that the problems of interest are in most cases underdetermined [12]. In fact, owing to the sophistication of numerical methods and increasing computational power, the forward model in diffuse optical tomography is computationally deployed using a (locally) linear model with a Jacobian matrix that has a number of columns that typically exceeds by far the number of rows, and this causes the associated inverse problem to be underdetermined. The main cause for this is the fact that the number of measurements is often far smaller than the number of nodes at which one wishes to reconstruct the optical properties. One thus balances increased resolution through finer mesh sizes against the numerical costs (e.g., round-off errors, computational time) of being underdetermined [13]. In this work the focus is on a particular numerical cost to being undetermined. While the emphasis here is on aspects relating to FT, the conclusions that are reached are also believed to be generally applicable to other applications such as DOT and bioluminescence tomography because of the similarity in their mathematical formulation [1416].
The results of this paper are as follows. It is shown, using standard structures in linear algebra, that the forward model in FT can project to its own kernel a significant portion of the vector containing the simulated optical properties in a manner that is not recoverable by standard regularization techniques such as Tikhonov regularization, truncated singular value decomposition, and methods based on Landweber–Fridman iterations. Here this portion of the optical properties vector (i.e., the component that is projected to the kernel) is referred to as “information loss,” and this loss is quantified using numerical experiments based on diffusion light transport modeling. It is important to state that this information loss is intrinsic to optical tomography with diffused light and that it is distinct from the propagation of data–model mismatch errors due to the ill-conditioned nature of the problem. In particular, it is shown here that the magnitude of the information loss can be as high as 50% of the magnitude of the optical properties vector for typical optical tomography geometries where a limited number of measurements are available [17]. The analysis presented here demonstrates that the loss of information increases with increasing node density using a finite element discretization scheme for the forward model based on the open source diffuse optical tomography simulation software package NIRFAST [18]. The main conclusion of the paper is that one must balance the desire for higher resolution against the accumulation of this measurable destruction of information. Though this loss of information is facilitated by the fact that the inverse problem in FT is typically underdetermined, it is demonstrated mathematically that being underdetermined alone is not the controlling factor. To this effect, we present evidence showing that the smoothing properties of the elliptic forward model that is used in the diffusion approximation to light transport in tissue may be partially responsible for the information loss.
The paper is organized in the following way. In Section 2 background material is presented, including a review of the diffusion equation used to model light transport (Subsection 2.A) and a brief overview of how the FT inverse problem is formulated as a linear matrix model (2.B). Then Subsection 2.C reviews the main linear algebra concepts required to understand the results presented. In Section 3 the methods and numerical experiments are described in detail, and Section 4 contains the results. The paper is concluded with discussion in Section 5 where an analysis of the results are presented as well as their implication for the field of tomography in highly scattering media. Appendix A contains some auxiliary facts from linear algebra.
A. Light Transport and Fluorescence Model
The forward problem for diffuse optical tomography may be described as follows: given a distribution of light sources on the boundary of a medium Ω and a distribution of optical parameters in the medium (absorption and scattering properties), find the resulting measurement set on the boundary. If the magnitude of the isotropic fluence is much greater than the anisotropic fluence, then it is well established that a diffusion equation model of light transport is an accurate reduction of the radiative transfer problem [5,6,19,20]. This diffusion approximation in the frequency-domain form is given by
equation M1
(1)
where Φ (r, ω)is the photon fluence rate, a position r [set membership] Ω and modulation frequency ω, S is an isotropic source term, equation M2 is the diffusion coefficient with absorption coefficient μa and reduced scattering coefficient equation M3. Finally, cm (r) is the speed of light in Ω at the point r defined by equation M4, where n(r) is the index of refraction at r and c is the speed of light in a vacuum. The air–tissue boundary is represented as a Robin boundary condition; that is, the flux leaving the boundary of the medium is the fluence rate at the boundary weighted by a factor that accounts for some internal reflection back into the tissue:
equation M5
(2)
where B is the term that deals with the refractive mismatch of the tissue with the surrounding media and ν is the unit surface normal of a patch on the boundary of the tissue.
The work in this paper is specifically focused on fluorescence molecular optical tomography in which fluorescent reporters are used to localize optical contrast via the accumulation or retention of the reporters; see, e.g., [8,21]. In this setting, a coupled pair of frequency-domain diffusion equations are used to predict the excitation and emission fluence rates Φx(r, ω) and Φm(r, ω), respectively. These equations are given as
equation M6
(3)
equation M7
(4)
where the diffusion coefficients are labeled x, excitation, and m, emission, and the fluorescence source term Sm(r, ω takes the form
equation M8
(5)
Here Qf is the quantum efficiency of the fluorophore, μa,f the absorption coefficient, and τ is the lifetime. Both the emission and excitation fluence equations satisfy a Robin boundary condition on [partial differential]Ω, as given in Eq. (2).
B. Fluorescence Tomography
In the limit where the absorption due to the fluorophores (μa,f) is much smaller than that of the other chromophores (e.g., hemoglobin, water, lipids), the solution to the coupled differential equations [Eqs. (3) and (4)] can be formally written as
equation M9
(6)
where G is the Green’s function associated with Eq. (4). Furthermore, Eq. (6) can be formulated on a discretized mesh of the volume Ω as
equation M10
where the integer k is used to label the spatial location of a measurement on the boundary [partial differential]Ω and the index i labels each of the Nν elements with a volume ΔVi in the mesh (e.g., nodes in a finite element volume).
Each fluorescence measurement in a full tomography data set is modeled using this discretized equation above. In this work, those solutions are obtained using the NIRFAST software [18] for which the interrogated domain is discretized using a finite element mesh. The fluorescence tomography problem is set up by considering several measurements around the periphery of the imaging domain. Those measurements and the discretized version of the corresponding differential equation are concatenated to form the matrix problem,
equation M11
(7)
where the measurement index k takes values up to the total number of measurements, Nm. The elements of the resulting Jacobian matrix, A = [aki], are
equation M12
(8)
In the remainder of this work, it is assumed that light transport physics in fluorescence tomography can be captured by this simple linear model, and we refer to the matrix A as either the reconstruction Jacobian or the forward model matrix, although they may differ in some cases where, for example, different meshes are used for forward modeling and image reconstruction.
C. Linear Algebra Primer: Singular Values and the Kernel Subspace
In this section the reader is reminded of certain standard linear algebra structures that we need to reference in order to establish the results in this paper. Of the multitude of possible references on the topic of linear algebra, we will use [22,23]. Only the basic definitions and results that are needed by us are established below; those who are looking for further details and arguments can find them either in the references or in Appendix A.
1. Vector Spaces and Linear Transformations
In all that follows, the symbols Nν [the number of nodes in a finite element method (FEM) mesh] and Nm (the number of measurements, i.e., source–detector pairs) are positive integers. We denote a standard real finite dimensional vector space by equation M13, where, for our purposes, k will take the values Nν or Nm. These vector spaces are endowed with the standard inner products [left angle bracket]·, ·[right angle bracket] given by [left angle bracket]x, y[right angle bracket] = xT y. The inner product is left unlabeled, and the context should inform the reader as to whether [left angle bracket]·, ·[right angle bracket] is an inner product on equation M14 or equation M15. In all instances in this paper, we will be working with a linear transformation L mapping vectors in equation M16 to vectors in equation M17; this is denoted by equation M18. The action of L can be represented by the action under matrix multiplication of a matrix A associated to L (see [22] for example); by abuse of notation we will use the symbol A for both a linear transformation and its matrix representation. Using standard notation, it is seen that equation M19, i.e., A has Nm rows and Nν columns. In the context of this paper, this matrix is the Jacobian matrix of the forward problem. Owing to the experimental design of diffuse optical tomography devices as well as the need for increased resolution, we will only consider the case where the dimension Nν of the domain vector space equation M20 is strictly greater than the dimension Nm of the range vector space equation M21; thus, we are explicitly dealing with an underdetermined system of linear equations.
The action of A splits equation M22 into the direct sum of two canonical subspaces. The kernel of A is a vector subspace of equation M23 consisting of all vectors that are mapped to the zero vector under the action of A; that is
equation M24
The kernel of A is always nonempty as 0 [set membership] Ker(A) for any choice of A. We let Ker(A)[perpendicular] be the orthogonal complement in equation M25 of Ker(A), i.e.,
equation M26
Note that Ker(A)[perpendicular] is itself a vector subspace of equation M27 . Then equation M28 is the direct sum of Ker(A) and Ker(A)[perpendicular], i.e.,
equation M29
(9)
This means that every vector equation M30 can be uniquely written as the sum of two vectors x0 and x[perpendicular], where x0 [set membership] Ker(A) and x[perpendicular] [set membership] Ker(A)[perpendicular]. It is clear that, if ν is a nonzero vector Ker(A)[perpendicular], then ≠ 0.
We recall that the range of equation M31 is the vector subspace of equation M32 defined by
equation M33
The rank of the linear transformation A, denoted by Rank(A), is the dimension of Ran(A); we remind the reader that the dimension of a vector space is the maximal number of linearly independent vectors in that space. The kernel and the range of A are tied together by the following fundamental theorem.
Theorem 2.1 (Rank Theorem [22])
Let A: equation M34 be a linear transformation. Then
equation M35
Given the standing assumption that Nν > Nm, and using the Rank Theorem, it is observed that the dimension of the kernel of A is always strictly positive since the range of A is contained in equation M36, and therefore rank(A) ≤ Nm. Finally, a useful piece of information concerning the action of A on the vector subspace Ker(A)[perpendicular] (see, e.g., [22]) is observed: the linear mapping A: Ker(A)[perpendicular] → Ran(A) is a bijection (i.e., a one-to-one and onto correspondence).
2. Singular Value Decomposition
There exists a change of coordinates, for both equation M37 and equation M38, so that the action of A in these coordinates is a diagonal action, that is, multiplication by a (nonnegative) constant in each coordinate direction. This change of coordinates, including the nonnegative constants, is called the singular value decomposition (SVD) of A (see, e.g., [23], Chapter 2). Under the assumption that Nm < Nν, the singular value decomposition of A = [aki], 1 ≤ kNm and 1 ≤ iNν, is given as f {νi, (uk), (σk)} so that (uk) is an orthonormal basis for equation M39,(νi) is an orthonormal basis for equation M40, and so that A(νk) = σkuk for 1 ≤ kNm with σk ≥ 0 and σk−1σk. Using matrix notation, we can represent equation M41 as A = UVT where U and V are Nm × Nm matrices constructed with columns from {uk} and {νi}, respectively, and ∑ is a Nm × Nν rectangular diagonal matrix with entries σj filling the diagonal.
The SVD provides a computationally convenient basis for the kernel of A (see Proposition 6.1 in Appendix A). It is with this basis that the magnitude can be computed of the projection by the FT forward model to its kernel.
A standard numerical approach to the study of the solutions to Eqs.s (3) and (4) is the FEM [18,24,25]; in fact, the FEM is used to approximate the forward problem concerning the fluence rate Φ in any given domain. Our FEM simulation is implemented via the software package NIRFAST [18,26], which is free-ware developed at Dartmouth for the purpose of simulating light propagation in biological tissue.
The examples considered here are based on simulations for two-dimensional triangular meshes of circular tissue-simulating phantoms of diameter 5 cm. The size of the phantom has been chosen to mimic small animals such as rats [27]. For the simulations, optical properties were chosen that were similar to what is expected in a small animal, namely equation M42 and μa = 0.03 mm−1. For simplicity, the optical properties were assumed to be the same at both excitation and emission wavelengths. The application of interest is the reconstruction of fluorescent sources representing tumors. As such, two inclusions of different sizes were simulated: one of radius 4 mm, and one of radius 2 mm, both with μa,f values ten times that of the simulated autofluorescence of the surrounding tissue. The images in the top row in Figs. 1 and and22 correspond to the fluorescence distribution (μa,f) in the phantom.
Fig. 1
Fig. 1
(Color online) Image reconstruction and information loss for an imaging geometry with 32 sources and 5 detectors. The top row shows the true solutions, the middle row shows the reconstructed images using the maximal number of SVD modes, and the bottom (more ...)
Fig. 2
Fig. 2
(Color online) Image reconstruction and information loss for an imaging geometry with 128 sources and 5 detectors. The top row shows the true solutions, the middle row shows the reconstructed images using the maximal number of SVD modes, and the bottom (more ...)
Imaging simulations were performed by varying two different types of parameters: detection geometry and FEM mesh resolution. First, different imaging geometries were considered that are consistent with the imaging system presented in [17,28]. Briefly, in those imaging systems, illumination is achieved with a pencil-beam laser source at different points on the surface. For each illumination point, focused detection is done in transmission using five photodetection channels. This detection geometry is installed on an imaging gantry that can rotate around the specimen allowing as many optical projections as desired to be acquired. In this work, numerical simulations were done associated with the following geometries: 8 sources (40 measurements), 16 sources (80 measurements), 32 sources (160 measurements), 64 sources (320 measurements), 128 sources (640 measurements), and 360 sources (1800 measurements). The other variable parameter that is considered in the simulations is the spacial resolution of the mesh. FEM meshes of the circular phantom were created associated to eight different degrees of spatial resolution corresponding to the following number of nodes: 106, 160, 275, 499, 736, 1207, 2333, 6414.
As explained in Section 2, diffuse fluorescence tomography can be treated as a linear problem in the limit that the absorption associated with the fluorophores is small relative to the absorption caused by tissue chromophores such as hemoglobin.
For each combination of FEM meshes and detection geometries, a Jacobian matrix A is calculated based on solutions to the diffusion equation; see Eq. (8). Then, using a basis for its kernel from the SVD, the projection of μa,f to the kernel is found. As explained below, this projection is directly associated with information loss that cannot be reconstructed using the standard least-squares regularization techniques presented in this work. This nonrecoverable portion of μa,f, i.e., the projection of μa,f to the kernel of the Jacobian, is displayed in the bottom row in Figs. 1 and and2,2, whereas the portion of μa,f that does lie in the orthogonal complement of the kernel (i.e., information that is in principle available for image reconstruction using the standard techniques described below in Section 4) is displayed in the middle row.
It is remarked that the portion of μa,f that lies in the orthogonal complement of the kernel (i.e., the middle row in Figs. 1 and and2)2) could also be obtained by computing a data vector ydata = a,f and then performing image reconstruction using the SVD decomposition based on the simulated data vector. The projection to the kernel of μa,f (i.e., the information loss) could then be computed by taking the difference between the original μa,f and the reconstruction. Note that this last statement applies only to the situation where no errors are introduced in the data vector prior to SVD reconstruction. For example, using different meshes for forward modeling and SVD image reconstruction would lead to another type of information loss due to the propagation of discretization and model mismatch errors in the reconstructed images.
In order to quantify the amount of information loss through projection to the kernel, we show in Fig. 3 the scaled error between the target image and the recovered image, computed using the equation
equation M43
as a function of mesh density. Here equation M44 denotes the reconstructed image, whereas equation M45 denotes the projection of μa,f to the kernel of A. The norm || · || is calculated as the standard Euclidean norm, that is, the square root of the sum of the squares of the values at all nodes. The mesh density (in units of nodes per millimeter) is calculated as the square root of the number of nodes divided by the area of the phantom.
Fig. 3
Fig. 3
(Color online) Reconstruction error versus node density for all detection geometries considered. The error decreases as the number of measurement increases, while the error typically increases with the number of nodes.
As has been proposed here, the displayed errors actually correlate with the size of the nonrecoverable information (as a proportion of the size of the true solution).
Our numerical experiments for those imaging geometries where the number of measurements is small compared to the number of nodes in the mesh show that the forward model A in FT projects to its kernel a significant portion of the optical properties vector, and this projection increases as the spatial resolution of the light transport mesh increases. This projection—regarded as a vector νo [set membership] Ker(A) in the canonical decomposition equation M46—has a large relative magnitude (often between 1/3 and 1/2 of the optical properties vector) and is unrecoverable via certain canonical regularization reconstruction methods. We therefore term this projection as “information loss.” Our numerical experiments using a diffuse light propagation model provide explicit examples demonstrating this loss of information in actual fluorescence reconstructions using an imaging geometry similar to what is commonly used in small-animal FMT [17]. We note that we are also considering, in our numerical experiments, detection geometries with a large enough number of source–detector pairs and a small enough number of nodes in the FEM mesh so as to no longer result in an underdetermined reconstruction problem. These correspond to those portions of the graphs in Fig. 3 that are flat and on the horizontal axis; i.e., there is no observed reconstruction error for these experiments. We have included these experiments in this exposition to demonstrate that, in order to increase the resolution in FT reconstruction, not only is a fine mesh necessary but in addition the detection geometry needs to be such that a large enough number of measurements exist to avoid information loss through projection into the kernel.
In what follows we examine mathematically the meaning of the projection to the kernel of the forward model and explain why this projection results in information loss. By analyzing certain standard regularization methods, we demonstrate that projection into the kernel of the FT forward operator always results in reconstructed images where the information contained in the kernel has been lost and thus cannot be recovered.
We begin with the following elementary example. A restricted form of the underdetermined least-squares inversion method can be formulated as follows: for a measurement vector y [set membership] Ran(A), the so-called least-squares solution μLS is given by
equation M47
that is, we choose for a solution μLS out of the collection of all possible solutions of = y as being the one that minimizes the square of the norm. Suppose now that a particular true optical properties distribution μtrue is known; that is, “true” has the meaning that, for a particular domain, the optical properties have been exactly measured. Under the forward model A the observation ymeas = true is calculated, and then, to test the efficacy of the least-squares inversion method under the assumption of this model A, a reconstruction is performed by calculating the least-squares solution μLS from the observation ymeas. Observe first that μtrue decomposes uniquely via the direct sum [Section 2.C, Eq. (9)] into equation M48, where μtrue,0 lies in the kernel of A and equation M49 is in the orthogonal complement of the kernel of A. Likewise, observe as well that any solution μ that satisfies = ymeas must be of the form μ = μ0 + μ[perpendicular], where μ[perpendicular] is in Ker(A)[perpendicular] and μ0 is in Ker(A). Note that, for any choice of solution μ in the set of all possible solutions, the projection of μ into Ker(A)[perpendicular] must be identical to equation M50, for otherwise A would not act in a one-to-one fashion on Ker(A)[perpendicular]. Using the fact that μ0 and equation M51 are orthogonal to each other for any choice of μ0 [set membership] Ker(A), then amongst all such solutions μ the minimum norm solution is equation M52. It is thus observed that the projection μtrue,0 of μtrue onto the kernel of A cannot be reconstructed in this setting, and in this sense μtrue,0 must be considered as being information that is “lost” by the forward model A under this form of the method of least squares.
In fact, for a given forward model the information associated with vectors contained in the kernel cannot always be retrieved given a particular choice of regularization. In support of this assertion, we consider several standard regularization techniques and demonstrate the phenomenon of information loss in optical tomography from first principles.
A. Tikhonov Regularization
More generally, let equation M53 be a linear transformation. The method of Tikhonov regularization (see [29,30]) seeks to mediate between fitting observed measurements equation M54 and fidelity to prior knowledge of some set of characteristics (e.g., size or smoothness) of the true solution equation M55. In its simplest form, Tikhonov regularization, for α [set membership] (0, ∞), is the solution μα of the operator equation (Theorem 2.5, [30]):
equation M56
(10)
which is derived as the solution to the minimization problem
equation M57
(11)
We work with Eq. (10):
equation M58
(12)
Recall that the range of the transpose AT is the orthogonal complement to the kernel of A (e.g., see [30], Appendix A). Thus, the right-hand side of Eq. (12) represents an element of the orthogonal complement of the kernel of A, implying that μα must lie in Ker(A)[perpendicular]. Thus, Tikhonov regularization, as realized in Eq. (11), cannot reconstruct that part of any solution that lies in the kernel.
Two more examples are now provided of popular reconstruction techniques that produce solutions in the orthogonal complement of the kernel of A.
B. Truncated Singular Value Decomposition
A brief overview of the SVD of the matrix A has been presented in Subsection 2.C.2. The SVD can also be used as a reconstruction technique. This method can be summarized as one that uses geometric information about the matrix A to project the measurement vector ymeas into the range of A. One then solves the equation = Proj(ymeas) in chosen a subspace of the orthogonal complement of the kernel of A and uses this value as the reconstructed solution of = ymeas.
A short description of this method is provided. Recall from our previous discussion of the SVD of equation M59 that there exist orthogonal matrices U and V and a diagonal matrix ∑, so that
equation M60
(13)
with
  • equation M61 and equation M62.
  • equation M63 is a matrix in diagonal form and whose diagonal elements satisfy
equation M64
The focus remains on the underdetermined case (Nm < Nν). Thus, the index value p, which denotes the smallest nonzero singular value equation M65, is strictly less than Nν; for ease of exposition we will assume that equation M66. Via the Rank Theorem (Theorem 2.1), the kernel of A is a nontrivial subspace of equation M67. In the general implementation of the truncated singular value inversion method, an index equation M68 is chosen that denotes the number of positive singular values to be used in the reconstruction; once again, for ease of exposition, we choose equation M69. Since the columns of U are a basis for equation M70, we have that equation M71. A consequence of Eq. (13) is that equation M72, equation M73 and equation M74, equation M75 (see Corollary 6.2 in Appendix A). Expressing both μ in terms of the basis {νj} and ymeas in terms of the basis {ui} and taking the inner product of both sides of = ymeas with ui, it is found that the truncated singular value solution equation M76 is
equation M77
which is in the orthogonal complement of the kernel of A since, as observed in Proposition 6.1 in Appendix A, {νp+1, …νNν} is a basis of the kernel of A.
C. Landweber–Fridman Iteration
Landweber–Fridman iteration is a reconstruction scheme that is based on fixed point iteration. The key point is that the iteration of a contracting map on a closed and complete space contains a unique fixed point and the computed fixed point is then taken to be the reconstruction. Using the normal equations
equation M78
(14)
on ymeas = , where P is the projection of ymeas onto the range of A in equation M79, one can show that the affine map
equation M80
is contracting for any choice of β so that equation M81, where σ1 is the largest singular value in the SVD of A. If we start at μ0 = 0, it is an easy observation that μk+1 = T(μK) is in Ker(A)[perpendicular] for all k. As Ker(A)[perpendicular] is a vector subspace of equation M82, the fixed point will be in Ker(A)[perpendicular].
The above reconstruction methods do not represent a complete list of all methods for which the projection to the kernel is unrecoverable; another example is Kaczmarz iteration as implemented in algebraic reconstruction technique [30].
It is well known that optical tomography images are plagued by poor spatial resolution. Past studies have linked this limitation to the physical nature of photon transport in tissue, i.e., that it is highly scattered. The highly ill-posed nature of the diffuse tomography models means that even minute levels of noise or model mismatch can push an inversion scheme away from finding biologically reasonable solutions. The choice and use of regularization techniques in optical tomography is thus a matter of import and delicacy. In this work new light has been shed on the optical tomography problem by investigating an intrinsic maximum information content, associated with a given optical tomography problem, that can be recovered using certain standard regularization techniques. We have demonstrated that this information content is intimately related to the nature of the imaging geometry and that it is also controlled by the spatial resolution of the mesh with which reconstruction is achieved. The analysis presented here has purposefully avoided the introduction of such issues as propagation of “errors” such as stochastic noise, general discretization error, and light transport model mismatch. The goal was to focus solely on a particular meshing error associated to the discretization of the FT forward model and to analyze numerically how it is affected by varying such parameters as the number of optical measurements and the spatial resolution of the FEM mesh used to generate the forward model. The sole optical imaging model studied in this work is diffuse fluorescence tomography, though we believe that our results are also applicable to the related diffuse optical tomography and bioluminescence tomography problems.
More precisely, we have been able to demonstrate, and quantify, that a significant portion of an optical properties vector is mapped by the FT forward model into its kernel, and examples have been given of regularization processes that cannot reconstruct the part of the vector that is mapped into the kernel. It has further been demonstrated numerically, under the assumption that the number of measurements stays fixed, that the size of the projection to the kernel grows as the mesh spatial resolution grows (number of nodes). In other words, the larger the dimension of the matrix kernel, the more information is potentially lost during the reconstruction process. It is true, as one increases the mesh density while leaving the number of measurements fixed, that the problem is becoming more underdetermined, but that in itself is not enough to mathematically guarantee that the size of the projection to the kernel will grow (see Example 6.3 in Appendix A). Thus, this property of projecting to the kernel in increasing amounts as the problem becomes more underdetermined is a fundamental property of the diffuse optical tomography forward model and not a mathematical artifact.
As mentioned earlier, we have found initial evidence that shows that the smoothing properties of the elliptic forward model that is used in the diffusion approximation to light transport in tissue may be partially responsible for the information loss. To corroborate this assertion, we show in Fig. 4 a graph that depicts the relative size of the projection to the kernel of optical properties vectors as a function of the smoothness of these vectors. The different smoothing states of the optical properties vectors were obtained by successively applying increasing amounts of Gaussian smoothing to a fixed initial vector. It is evident that the smoother the optical properties vector is, the smaller the projection into the kernel of the forward operator becomes. We plan to explore this phenomenon further in a subsequent paper.
Fig. 4
Fig. 4
(Color online) Reconstruction error as a function of smoothness of the optical properties vector. A value of 1 indicates that no smoothing has been applied (and thus the optical properties vector exhibits sharp edges), whereas a value of 10 indicates (more ...)
The numerical experiments presented in Fig. 3 also demonstrate that the magnitude of the information loss (size of the projection to the kernel relative to the size of the projection to its complement) depends on the detection geometry that is considered, with preferable geometries being associated with a larger number of optical measurements. An interesting case is the 8 sources (40 measurements) curve for which the reconstruction error remains relatively flat for all node densities. In this case, because there are so few measurements, the reconstruction error remains high even when the number of nodes is at its smallest. It should be noted in closing that the work presented here is expected to have a significant impact mainly for tissue imaging geometries where the number of measurements is limited, such as point-by-point detection configurations based on photomultiplier tubes. In fact, the impact of information loss through the forward matrix kernel is not expected to be as significant when high density of measurements are available using multiple projections based on wide-field CCD detection. Several publications have demonstrated in the past that using high densities of measurements can afford significant spatial resolution gains for both FT [31] and DOT [32] instruments.
An important conclusion that comes out of the analysis presented here is that, in situations when the number of measurements is limited, one must balance the search for increased spatial resolution against the increased loss of information by the FT forward model under an increasingly fine meshing of the domain. For a fixed number of sources and detectors, as one increases mesh density in an attempt to improve resolution, it is important that one avoids certain standard regularization techniques that, for mathematical reasons, cannot reconstruct the information that is lost by the forward model. In fact, it was demonstrated here that, if image reconstruction in optical tomography is to be done using regularization techniques such as Tikhonov, truncated SVD, or Landweber–Fridman, for example, then in order to gain a finer resolution it is not sufficient to increase the density of the reconstruction mesh, and the only available option is to increase the number of measurements. As a rule of thumb, meshes for which the number of nodes is comparable to the number of measurements will suffer less from the information loss described in this paper, and typically no gains can be expected from using significantly higher-resolution meshes. Of course, there are situations where high-resolution meshes are unavoidable when, for example, sharp boundaries need to be modeled. However, in such situations if possible the number of sources and detectors should be increased to offset the information loss associated with the higher resolution. It is of course well known that increasing the number of measurements is beneficial in optical tomography (see, e.g., [32,33]). The findings here are rather that, under limited tissue sampling conditions, increasing the spatial resolution of the reconstruction mesh can be detrimental to image quality and that this is not only caused by the underdetermined nature of the problem. Further and preliminary numerical experimentation has also revealed that boundary regions in a true solution (e.g., edges of a fluorescent lesion) account for a significant portion of the projection to the kernel. We hope to explore this issue, along with possible inverse technique remedies, in a subsequent paper.
ACKNOWLEDGMENTS
This work has been funded by National Science Foundation grants DMS-1031954 and DMS-1031955 and National Cancer Institute grants RO1 CA120368 and K25 CA138578. The authors P. B.-T. and E. C. T. thank the Thayer School of Engineering at Dartmouth College for providing a stimulating and conducive environment for research during the academic year 2010–2011.
APPENDIX A. SINGULAR VALUE ANALYSIS AND AN EXAMPLE
This appendix has two purposes. First, in service of making the paper self contained, a pair of standard linear algebra results used in the paper is established. Second, an elementary example is provided of a matrix having a high-dimensional kernel along with a vector of large magnitude that does not project into the kernel. The point of this observation is that being underdetermined is in and of itself not enough to guarantee information loss; rather it is a property of the FT forward model itself.
The reader is reminded of the standing assumption that Nν > Nm, and we let equation M83 be the smallest positive singular value. Both results that are established are independent of the condition NmNν, but both to avoid switching notation and because we are only concerned with the undetermined case, this convention is maintained. The following proposition gives a computationally convenient basis for the kernel of A, in terms of the SVD, that we used to calculate the size of the projection into the kernel.
Proposition 6.1
Let equation M84. Then the collection of orthonormal vectors {νp+1, ….., νn}, chosen from (νi), 1 ≤ iNν, is an orthonormal basis for the kernel of A.
Proof
Owing to the assumptions that σp is the smallest non-zero singular value and that the collection (uj) is orthonormal, we can assert that the range of A contains a maximal collection of p linearly independent vectors, which implies that the rank of A is p. Thus, via the Rank Theorem, we have that the dimension of the kernel of A is equation M85.
Let equation M86. Then equation M87. Since the collection (νi) (1 ≤ i) is orthonormal, we know that equation M88 for all ik and furthermore that equation M89 for all 1 ≤ kn. The columns of the matrix V consist of the vectors vi, and thus VT vk is a vector all of whose components are zero, except for the kth component, which equals 1. But multiplying this vector from the right with the diagonal matrix ∑, whose diagonal entries from place p+1 on (and in particular in place k) are zero, results in the zero vector. This implies that Avk = 0, so that each orthonormal vector in the collection {νp+1, …, νNv} is in the kernel of A. As any collection of equation M90 linearly independent vectors in a equation M91 dimensional vector space is a basis, we observe that {νp+1, …, νNv} is a basis of the kernel of A.
A useful and well-known corollary to this proposition is recorded. The SVD of the matrix A also provides a change of basis that reveals that the transpose AT also has a diagonal action.
Corollary 6.2
Let {(νi), (uj), (σj)}, be the SVD of A. Then for 1 ≤ jp, we have
equation M92
Proof
equation M93 of which {ν1, ….., νNv} is a basis, and so
equation M94
Finally, we finish with the promised example. Fix a positive integer Nm and let {Nvi} be a sequence of integers greater than Nm for which limi→∞Nνi = ∞. Furthermore, let equation M95. Assume as well that Ai is of rank Nm for all i. The Rank Theorem would thus tell us that, for each i, the dimension of the kernel of Ai is NνiNm (and thus going to infinity as i increases). Recall that Nνi represents the number of nodes in our mesh at which we wish to reconstruct optical properties, whereas Nm is the number of measurements (observations) that we have made. We have shown in this paper, for in the FT tomography setting, that as Nvi increases, the size of the projection onto the kernel of our optical properties (and therefore the model inherent information loss) also increases.
This is an intrinsic property of the FT forward model and not a mathematical artifact, as the following simple example will show.
Example 6.3
The matrix Ai (i ≥ 2), as a linear transformation, is in equation M96 and explicitly has the form that the first row has ones in the odd numbered columns with zeroes in the even numbered columns while the second row has zeroes in the even numbered columns and ones in the even numbered columns. That is,
equation M97
Once again it is easy to see that the rank of Ai is two, regardless of the choice of i, and so by the Rank Theorem the dimension of the kernel of Ani is 2i – 2 for each choice of i. In this example it is easier to focus on a particular example, as the general case follows easily from induction arguments. Thus, we will explicitly consider the i = 3 case A6, i.e.,
equation M98
Explicitly calculating, we find a basis for the kernel of A6 as
equation M99
where here the idea is to first show that the collection of vectors is linearly independent and then note, for any finite dimensional vector space of dimension k, that any collection of k linearly independent vectors is indeed a basis.
It is now an easy observation that the constant vector x = (1, 1, 1, 1, 1, 1) is orthogonal to each of the basis vectors, and thus the projection of x onto the kernel of A6 is the zero vector (that is, “no information” in (1, 1, 1, 1, 1, 1) is lost to the kernel by the action of A6).
Footnotes
OCIS codes: 110.6960, 110.7050, 110.3010, 100.6950.
1. Shah N, Cerussi A, Eker C, Espinoza J, Butler J, Fishkin J, Hornung R, Tromberg B. Noninvasive functional optical spectroscopy of human breast tissue. Proc. Natl. Acad. Sci. USA. 2001;98:4420–4425. [PubMed]
2. Pogue BW, Poplack SP, McBride TO, Wells WA, Osterman KS, Osterberg UL, Paulsen KD. Quantitative hemoglobin tomography with diffuse near-infrared spectro-scopy: pilot results in the breast. Radiology. 2001;218:261–266. [PubMed]
3. Zhang Q, Brukilacchio TJ, Li A, Stott JJ, Chaves T, Hillman E, Wu T, Chorlton M, Rafferty E, Moore RH, Kopans DB, Boas DA. Coregistered tomographic x-ray and optical breast imaging: initial results. J. Biomed. Opt. 2005;10:024033. [PubMed]
4. Arridge SR, Hebden JC. Optical imaging in medicine: II. Modelling and reconstruction. Phys. Med. Biol. 1997;42:841–853. [PubMed]
5. Arridge SR. Optical tomography in medical imaging. Inverse Probl. 1999;15:R41–R93.
6. Jacques SL, Pogue BW. Tutorial on diffuse light transport. J. Biomed. Opt. 2008;13:041302. [PubMed]
7. Ntziachristos V, Tung CH, Bremer C, Weissleder R. Fluorescence molecular tomography resolves protease activity in vivo. Nat. Med. 2002;8:757–761. [PubMed]
8. Leblond F, Davis S, Vald’s P, Pogue B. Pre-clinical whole-body fluorescence imaging: review of instruments methods and applications. J. Photochem. Photobiol. B. 2010;98:77–94. [PubMed]
9. Leblond F, Tichauer KM, Pogue BW. Singular value decomposition metrics show limitations of detector design in diffuse optical tomography. Biomed. Opt. Express. 2010;1:1514–1531. [PMC free article] [PubMed]
10. Leblond F, Dehghani H, Kepshire D, Pogue BW. Early-photon fluorescence tomography: spatial resolution improvements and noise stability considerations. J. Opt. Soc. Am. A. 2009;26:1444–1457. [PubMed]
11. Culver JP, Ntziachristos V, Holboke MJ, Yodh AG. Optimization of optode arrangements for diffuse optical tomography: a singular-value analysis. Opt. Lett. 2001;26:701–703. [PubMed]
12. Arridge SR, Lionheart WRB. Nonuniqueness in diffusion-based optical tomography. Opt. Lett. 1998;23:882–884. [PubMed]
13. Yalavarthy P, Dehghani H, Pogue B, Paulsen K. Critical computational aspects of near infrared circular tomographic imaging: analysis of measurement number, mesh resolution and reconstruction basis. Opt. Express. 2006;14:6113–6127. [PubMed]
14. Contag CH, Ross BD. It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology. J. Magn. Reson. Imaging. 2002;16:378–387. [PubMed]
15. Contag CH, Bachmann MH. Advances in in vivo bioluminescence imaging of gene expression. Annu. Rev. Biomed. Eng. 2002;4:235–260. [PubMed]
16. Dehghani H, Davis SC, Jiang S, Pogue BW, Paulsen KD, Patterson MS. Spectrally resolved bioluminescence optical tomography. Opt. Lett. 2006;31:365–367. [PubMed]
17. Kepshire DS, Mincu N, Hutchins M, Guber J, Dehghani H, Hypnarowski J, Leblond F, Khayat M, Pogue BW. A microcomputed tomography guided fluorescence tomography system for small animal molecular imaging. Rev. Sci. Instrum. 2009;80:043701. doi:10.1063/1.3109903. [PubMed]
18. Dehghani H, Eames M, Yalavarthy P, Davis S, Srinivasan S, Carpenter C, Pogue B, Paulsen K. Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction. Commun. Numer. Meth. Eng. 2009;25:711–732. [PMC free article] [PubMed]
19. Groenhuis RAJ, Ferwada HA, Ten Bosch JJ. Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory. Appl. Opt. 1983;22:2456–2462. [PubMed]
20. Groenhuis RAJ, Ten Bosch JJ, Ferwada HA. Scattering and absorption of turbid materials determined from reflection measurements. 2: Measuring method and calibration. Appl. Opt. 1983;22:2463–2467. [PubMed]
21. Paithankar D, Chen A, Pogue B, Patterson M, Sevick-Muraca E. Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media. Appl. Opt. 1997;36:2260–2272. [PubMed]
22. Friedberg S, Insel A, Spence L. Linear Algebra. 4th ed Prentice Hall; 2003.
23. Golub G, van Loan CF. Matrix Calculations. Johns Hopkins University; 1989.
24. Paulsen K, Jiang H. Spatially varying optical property reconstruction using a finite element diffusion equation approximation. Med. Phys. 1995;22:691–701. [PubMed]
25. Arridge SR, Schweiger M, Hiraoka M, Delpy DT. A finite element approach for modeling photon transport in tissue. Med. Phys. 1993;20:299–309. [PubMed]
26. Zhu Q, Dehghani H, Tichauer KM, Holt RW, Vishwanath K, Leblond F, Pogue BW. A three-dimensional finite element model and image reconstruction algorithm for time-domain fluorescence imaging in highly scattering media. Phys. Med. Biol. 2011;56:7419–7434. [PubMed]
27. Leblond F, Tichauer KM, Holt RW, El-Ghussein F, Pogue BW. Toward whole-body optical imaging of rats using single-photon counting fluorescence tomography. Opt. Lett. 2011;36:3723–3725. [PubMed]
28. Tichauer KM, Holt RW, El-Gussein F, Zhu Q, Dehghani H, Leblond F, Pogue BW. Imaging workflow and optical data calibration for CT-guided whole-body time-domain fluorescence tomography. Biomed. Opt. Express. 2011;2:3021–3036. [PMC free article] [PubMed]
29. Hansen PC. Discrete Inverse Problems: Insight and Algorithms. SIAM; 2010.
30. Kaipio J, Somersalo E. Statistical and Computational Inverse Problems. Springer-Verlag; 2005.
31. Chen J, Venugopal V, Lesage F, Intes X. Time-resolved diffuse optical tomography with patterned-light illumination and detection. Opt. Lett. 2010;35:2121–2123. [PubMed]
32. Konecky SD, Panasyuk GY, Lee K, Markel V, Yodh AG, Schotland JC. Imaging complex structures with diffuse light. Opt. Express. 2008;16:5048–5060. [PMC free article] [PubMed]
33. Wang ZM, Panasyuk GY, Markel VA, Schotland JC. Experimental demonstration of an analytic method for image reconstruction in optical diffusion tomography with large data sets. Opt. Lett. 2005;30:3338–3340. [PubMed]