1. Introduction

Conventional imaging modalities such as magnetic resonance imaging (MRI) and x-ray computer-aided tomography (x-ray CT) provide high-resolution medical imaging enabled by the direct geometric correlation between incident and detected radiation. Yet the high cost of operating hospital MRI facilities, and the inability to detect important diseases through x-ray CT imaging suggests an opportunity for the development of near-infrared (NIR) biomedical optical imaging. NIR biomedical optical imaging, or optical tomography, depends upon the low absorbance, yet high scattering of non-mutagenic near-infrared light in the “therapeutic wavelength window” (600–1000 nm) enabling it to safely propagate through several centimeters of tissues. While the propagation of low-energy NIR light can occur safely, the high scattering properties of tissues renders the direct geometric correlation between incident and detected irradiation invalid. The forward imaging problem, (i.e., prediction of the propagation of NIR light through tissues when a map of the tissue optical properties is known) can be described in terms of the diffusion approximation to radiative transfer. The inverse imaging problem (i.e., prediction of the interior optical properties from measurements of light propagation made at the exterior tissue-air interface positions) requires solution of a series of equations that are non-linear in the optical properties to be estimated. In this contribution, we describe a novel algorithm to estimate a solution to the nonlinear, inverse imaging problem for fluorescence frequency-domain photon migration (FDPM). In the following sections of this introduction, we briefly (i) present the background of frequency-domain photon migration (FDPM) imaging, (ii) compare prior work and our current approach towards solution of its inverse imaging problem, and (iii) introduce the concept of fluorescent contrast agents which can accelerate the convergence of the inverse imaging solution.

1.1 Frequency-domain photon migration: forward and inverse problems

Frequency-domain photon migration, depends upon launching intensity modulated (30-200 MHz) light at the interface of a highly scattering medium (such as tissue), and detecting the intensity-modulated wave that successfully propagates to the detector located a distance ρ away from the incident source. Depending upon the spatial distribution of interior absorption and scattering optical properties, the detected light is both phase-delayed by θ, and amplitude attenuated by factor M when compared to the incident light (see Figure 1). The solution to the forward imaging problem involves predicting the phase-delay, θ, and amplitude demodulation, M as a function of angular modulation frequency, ω, and position, r, along the tissue surface given a known spatial distribution of absorption and scattering optical properties. The absorption coefficient, μ_ax , representing the inverse mean absorption path and the isotropic scattering coefficient, μ′_sx, equivalent to the inverse isotopic scattering length, govern the propagation of light through scattering media, such as tissues. Since near-infrared light is multiply scattered in tissues, the numerical solution to the diffusion approximation of the radiative transfer equation provides solution to the forward biomedical optical imaging problem. For a propagating intensity modulated wave of light, the optical diffusion equation is written as

where Φ_x represents the complex number describing the scalar flux of photons at position r; ω is the modulation frequency; c is the speed of light within the medium; μ_axi is the absorption due to chromophores and μ_axf is the absorption due to fluorophores (μ_ax = μa_xi + μ_axf). D_x is an optical diffusion coefficient equivalent to 1/[3(μ_ax + μ′_sx)]. From numerical solution of Eqn (1) with a Robin boundary condition, the complex scalar flux at the surface, Φ_x = (=Me ^iθ) and the phase delay, θ, and amplitude attenuation, M, can be determined. It is important to note that the diffusion equation applies when μ_ax ≪ μ′_sx and D_x is approximated by 1/[3(μ′_sx)].

 

Figure 1 Schematic of fluorescence photon migration.

Download Full Size | PDF

The successful implementation of biomedical optical imaging involves the effective solution of the inverse imaging problem, i.e., determining the interior optical property map of absorption and scattering given measurements of θ and M along the tissue surface. To date, approaches have focused upon the linearization of the problem using first order Born or Rytov approximations [1, 2]. Non-linear optimization employing iterative Born [3, 4], or distorted Born methods [3,5] require accurate information regarding the normal “background” optical properties which the tissue “heterogeneity” resides. Since the normal tissue background optical properties can be expected to vary greatly, this a priori information is not realistically feasible for biomedical optical imaging.

Other investigators have used the modified Newton-Raphson method with the Levenberg-Marquardt regularization as the central inversion step [6, 7, 8 and 9]. While not requiring accurate “background” optical properties for initial parameter estimate guesses for successful reconstruction, the Newton-Raphson with Levenberg-Marquardt regularization has three distinct disadvantages for biomedical optical imaging. The first disadvantage is that the full Hessian matrix is not fully considered making the method inefficient for highly nonlinear problems such as optical imaging. Numerical experiments by McKeown [10] and Veapucci [111] have shown that while the Gauss-Newton and Levenberg-Marquardt methods solve nonlinear least squares problems efficiently when the residuals are small or zero at the solution, these methods are less efficient than the quasi-Newton method when the residuals are significantly nonlinear. Meyer [12] provides theoretical support of this observation. In the work presented herein, we demonstrate the use of a truncated Newton method with trust region for the efficient solution of the large scale, non-linear optical imaging problem. The second disadvantage to Newton-Raphson approaches is that the methods require expensive numerical calculation of the Jacobian that can be subject to error. In this work, we adapt modifications of the reverse automatic differentiation to compute the gradient defining the search direction efficiently and without numerical error [13, 14]. Finally, the Levenberg-Marquardt method requires storage of the full Jacobian matrix, which represents a limitation for optical imaging in three dimensions and when the number of unknowns increases. Since the truncated Newton method requires only the product of the Hessian with the direction at any point, the method is comparatively more efficient for large problems. Our inversion strategy is presented in its entirety in Section 3.

While the difficulties associated with the nonlinear and large scale of inverse imaging problem can be addressed with our inversion strategy, additional significant challenges exist when attempting to reconstruct absorption and scattering properties that are physiologically feasible. In a study involving actual optical property measurements of normal and diseased tissues, Troy et al. [15] found that the optical properties of absorption and isotropic scattering of normal and diseased breast tissues, for example, may not provide sufficient contrast for successful image reconstruction. Indeed, the small influence of “heterogeneity” absorption and scattering on FDPM measurements measured at the boundary can often be within or close to the measurement error of 1 degree in phase and 1% in amplitude attenuation.

1.2 Fluorescence frequency-domain photon migration: forward and inverse problems

In order to over come the difficulties associated with insufficient contrast and costly nonlinear optimization, we and other investigators have developed fluorescence FDPM imaging approaches and inversion strategies [8, 16, 17, 18]. Fluorescence FDPM physically depends upon the administration of a fluorescent contrast agent, or another agent or gene vector, which results in the expression of a fluorescent protein that emits in the near-infrared wavelength regime. When activated by the absorption of an intensity modulated wave, Φ _x (r, ω), an intensity modulated fluorescent wave, Φ _m (r, ω), is generated at the same modulation frequency, but is amplitude attenuated and phase-delayed relative to the activating excitation wave. The phase-delay, θ_m, and amplitude attenuation, M_m, of the generated fluorescent wave can be as high as 90 degrees and a factor of 10’s –100’s of the incident excitation wave owing to the fluorescence properties (lifetime and quantum efficiency) of the dye. Specifically, as the biochemical environment of the fluorophore changes, its lifetime, τ (i.e., the lifetime of the activated fluorophore, or the mean time between absorption of excitation and release of fluorescent light) also alters, providing discrimination of diseased tissues possibly on the basis of biochemical environment. The detected intensity modulated fluorescence wave that has been generated within and which has propagated to the tissue boundary possesses a phase-delay, θ_m, and amplitude attenuation, M_m, that is exquisitely more sensitive to embedded “heterogeneities” than possible with absorption FDPM imaging described above [16].

Predictions of fluorescence FDPM measurements of phase-delay, θ_m, and amplitude attenuation, M_m are achieved through the solution of the complex fluorescence fluence, Φ _m (r, ω), at the boundary from the fluorescence optical diffusion equation:

where D_m is the optical diffusion equation at the emission wavelength (=1/3⌊μ_am+μ′_sm⌋)≈1/3⌊μ′_sm⌋; μ_am and μ′_sm are the absorption and isotropic scattering coefficients at the fluorescent wavelength; and the right hand term describes the generation of fluorescence within the medium. The term ϕ represents the quantum efficiency of the fluorescence process and the absorption owing to fluorophore is represented by the coefficient μ_ax→m. Here we consider first order relaxation decays only. Note that the source term requires coupling with the solution of excitation fluence described by Eqn (1), where by the absorption coefficient at the excitation light, μ_ax is now provided by the naturally occurring, non-fluorescent chromophores, μ_axi, and the fluorescent contrast agent, μ_ax→m. The numerical solution for the excitation fluence distribution Eqn (1) with the Robin boundary condition enables prediction of Φ _m from Eqn (2) at the medium boundary and determination of θ_m and M_m. It is important to note that the diffusion equation applies when μ_am<<μ′_sm.

With the added unknowns at the emission wavelength, the disadvantages using Newton-Raphson optimization for fluorescence FDPM imaging may become more severe. Paithankar et al., [8] and Troy and Sevick-Muraca [19] employed Newton-Raphson, multi-grid finite difference reconstructions of μ_ax→m τ, and ϕμ_ax→m, given correct priors on other optical properties. Recently, Jiang [20] performed similar work, employing dual meshing finite element methods instead of finite difference methods. Both Paithankar et al. [8] and Jiang [20] resorted to removal of “spurious” optical property values or “filtering” using empirically chosen parameters to achieve the “correct” image distorted through inefficient numerical computation of the Jacobian. O’Leary et al., [17] and Chang, et al., [18] employed Born and iterative Born which again not only required known “background” optical properties, but which also was unable to handle fluorescence in the “background.” Later Chang et al., [21] demonstrated the ability to reconstruct fluorescence lifetime under conditions of “imperfect” uptake given the spatial distribution of absorption owing to fluorophore; and Lee and Sevick-Muraca [22] employed an iterative Born-type solution of simultaneous absorption and lifetime from fluorescence FDPM synthetic data.

1.3 Fluorescence FDPM imaging: current work

In this contribution, we also focus upon fluorescence FDPM imaging, owing to the significant influence of interior fluorescent heterogeneities on measured signals at the medium boundary and adapt our inversion strategy for the non-linear, large scale optimization required for biomedical optical imaging of absorption and fluorescence lifetime. In the following section, the numerical methodology for solving the FDPM forward imaging problem and the simulator used for generation of synthetic data sets for testing our inversion algorithm are presented. Section 3 describes our formulation of the inverse imaging problem as an optimization problem using truncated Newton’s methods with trust region and the reverse automatic differentiation in order to reduce the computational overhead associated with optimization. The performance of our algorithm using synthetic data sets is presented in Section 4 and the conclusions, prognosis for FDPM imaging and ongoing work are summarized in Section 5.

2.0 Model prediction of FDPM

2.1 Boundary conditions for FDPM forward simulator

Prediction of FDPM measurements on surface dΩ resulting from diffusion of intensity-modulated waves within volume Ω, is accomplished from solution of Eqn (1) and (2) using the Robin boundary condition. We employ the Robin boundary condition of the form given by Ishimaru[23, 24]:

Where S is the complex flux density associated with a point source of excitation light located at r_s on dΩ; n is the unit vector normal to the surface, and γ accounts for refractive index mismatch at the boundary which through Snell’s reflection, results in internal reflection at the surface.

Groenhuis et al. [25] have developed an empirical approach to determine the value of the parameter γ, viz.

where

n_rel is the relative index of refraction of the scattering medium with respect to the surrounding medium.

2.2 Galerkin finite element formulation for solution of Φ_x

We employ the Robin boundary condition in our numerical solution of Φ_x,m employing the Galerkin finite element model with linear triangular elements. In our formulation, we assume μ_ax,m<<μ_sx,m, so that D_x,m ≈1/3μ_sx,m. The formulation of the Galerkin finite element solution of Φ_x begins by multiplying Eqn (1) by a weighting function, w_j, and integrating over the domain of interest:

where w_j {w_j : j = 1,⋯,N} is a set of linearly independent weighting functions. N is the number of unknows.

Using Green’s theorem, equation (6) becomes

where the second term is a line integral along the external boundary and w_j is assumed to be continuous over the whole domain. Now introducing equation (3) in the line integral

Upon combing Eqns (8) and (7), once can write

Suppose that the solution domain Ω is divided into M ₁ triangular elements. Now we can obtain the element equation from Eqn (9) as follows:

Assuming that ${\mathrm{\Phi }}_x^{\mathit{\text{el}}}$ , ${\mathrm{D}}_x^{\mathit{\text{el}}}$ and μ^el_axf vary linearly within each triangular element el and μ_axi is constant, we may then write ${\mathrm{\Phi }}_x^{\mathit{\text{el}}}$ , ${\mathrm{D}}_x^{\mathit{\text{el}}}$ and μ^el_axf for each element as:

where the L_i are the natural co-ordinates for the triangle [26]. According to the Bubnov-Galerkin method, the weighting functions are chosen to be the same as the approximation functions used to represent ${\mathrm{\Phi }}_x^{\mathit{\text{el}}}$ , that is,

The equation (10) can be written as:

After integration of Eqn (12) and rearrangement, the resulting element stiffness equations become:

Assembly of these element stiffness matrices ${\mathbf{K}}_i^{\mathit{\text{el}}}$ and local vector r ^el yields the system equations KΦ̄ _x =b [26].

Since S denotes the complex point source applied to the surface dΩ. at position r_s and S is expressed as:

where r_s is the position vector of the point source. By definition, we have

where el is the element of the finite element mesh. From Eqn(12) we have

If the position r_s of the point source coincides with the position of node j, when L_j(r_s) = 1, i.e. we have

2.3 Galerkin formulation of finite element solution Φ_m

The forward finite element solution of the emission diffusion equation is formulated similarly to the excitation diffusion equation with the exception that (i) the complex excitation fluence must be first computed, and that (ii) the formulation of fluorescence lifetime poses some additional considerations. As in the formulation of the solution for Φ_x, parameters of μ_ax→m, τ, and D_m are expected to vary linearly within the triangular elements. In order to facilitate the formulation of the Galerkin finite element method, the RHS term of Eqn (2) is expressed in terms of a binomial expansion assuming ωτ<1 and higher order terms are neglected. Hence, the RHS source term of Eqn (2) becomes S_m (r,ω)=ϕμ_ax→m[1 + iω(r)τ]Φ _x (r, ω). This assumption is appropriate, since the phase-delay, θ_m, of the generated fluorescent intensity modulated wave reaches a maximum of 90 degrees with respect the local excitation intensity modulated wave and becomes insensitive to fluorescence lifetime when ωτ>1. The formulation of the emission finite element solution for Φ _m is identical to that presented in section 2.2 for Φ_x with the exception that the parameters Φ _m , D_m, μ_ax→m, and τ are assumed to vary linearly within each element. The finite element formulation of emission diffusion is given below:

It is important to emphasize that the absorption owing to fluorophore is equivalent in our case to μ_ax→m and to μ_axf. We differentiate between the two parameters when reconstructing absorption from synthetic excitation, μ_axf, and absorption from emission data μ_{ax → m}.

2.4 Solution of the system of finite element equations

For both excitation and emission finite element formulations, the local stiffness matrices, ${\mathbf{K}}_{\mathrm{j}}^{\mathit{\text{el}}}$, are formed and then combined into a complex, global stiffness matrix, K, which depends upon the absorption and scattering properties at each wavelength. For solution of excitation fluence, Φ̄_x, the vector b reflects the complex source of excitation light from boundary point sources at positions r_s, while for solution of the emission fluence, Φ̄_m, vector b represents the source generated from within volume Ω.

To solve for Φ̄_x,m , the complex set of sparse, linear equations generated by the finite element method,

are efficiently solved using a subroutines ZGBTRF and ZGBTRS (LAPACK subroutines). ZGBTRF computes a LU factorization of a complex M ₂ by N band matrix using partial pivoting with row interchanges. (M ₂ = 2*KL + KU + 1, KL and KU are the numbers of subdiagonals and superdiagonals within the band of K, respectively). ZGBTRS solves a system of linear equations with a general band matrix K using the LU factorization computed by ZGBTRF. The main advantage of this method is that once the decomposition is done, the solution vector can be obtained for any number of right hand side vectors. The accuracy of the solution for Φ̄_x,m depends crucially on the mesh size, h, with the accuracy improving with h ².

3.0 Inverse solution for fluorescence FDPM

The variables in the inverse problems are the absorption coefficients, μ_axf (for excitation), μ_ax→m (for emission), and the lifetime τ at each nodal point of a finite element model. If there are N nodal points then the number of measurements must be greater than N. The number of measurements made at N_B boundary nodes for N_s sources, requires that

Due to the reciprocality of sources and detectors, i.e., when a detector position is swapped with a source position equivalent measurements are registered, then we require:

If l =1,…,N_s denotes the different distribution of source s and j =1,…, N_B denotes the boundary nodes at which measurements are made, the error function for optimization of absorption and lifetime imaging may be taken as

The subscript c denotes the values calculated by the forward simulator problem and the subscript me denotes the experimental (or in this case synthetically generated) fluence values. The superscript * denotes the complex conjugate of the complex function Φ_x,m. The subscript l and j denote fluence values that arise between source l and detector j.

It is seen from the equation (20) that the calculation of error function E involves N_s solutions of the direct problem. The global stiffness matrix K _x,m for excitation and emission remains constant with only the vector b changing with each excitation source. Since the global stiffness matrix K _x depends on μ̄_axf and b is independent of μ̄_axf, LU decomposition of K _x for solution of Φ_x need only be performed once in each iteration of the optimization problem.

An efficient way of calculating E_x,m for excitation and emission equations consists of the following steps:

As most of the efficient optimization algorithms, the truncated Newton method requires the function value E_x,m and the gradient of the error function, ∇E_x,m , in order to efficiently update parameter estimates. The gradients of the error function in Eqn (20) represent the direction in which the error function is increasing/decreasing most rapidly. Since absorption coefficients μ_axf μ_ax→m and τ are the unknown parameters, the gradients are obtained by taking partial derivatives of Eqn (20) with respect to μ_axf, μ_ax→m and τ.

〈.,.〉 denotes the inner product of two complex vectors and Re(.) denotes the real part of a complex. number. The gradient of the real valued error function with respect to real valued μ̄_axf μ_ax→m , and τ remains real valued. At minimum these gradients should be zero, so the real part of the gradients are zero. Therefore, we consider only the real part of the gradients. In the following, the optimization strategy using the truncated Newton method with the gradients computed by reverse automatic differentiation is described.

3.1 Truncated Newton’s method as an optimization strategy for FDPM imaging

The truncated Newton method with trust region is used for non-linear optimization problem. The reader interested in background literature is referred to the references [27, 28].

Newton’s method is based upon approximating the function E_x,m (μ̄_ak +d) (μ_a denotes μ_axf {for excitation} or μ_ax→m {for emission} in this section) at the kth iteration by the quadratic model:

In these equations g _k =g(μ̄_ak =∇E_x,m (μ̄_ak and G _k =∇² E_x,m (μ̄_ak) denote the gradient vector and Hessian matrix, respectively. The Newton direction is obtained from an exact minimum of above equation; i.e. the search direction d at iteration k is defined by the equation

A sequence of approximate solutions of the above equation is generated by the conjugate gradient method until a vector d _i (see Step 3) is obtained for which the following condition on the relative residual is satisfied as suggested by Dembo and Steihaug [29];

where r _i = G _k d _i + g _k. The conjugate gradient iteration is then truncated and d _i is used as the search direction for the minimization. The relative residual is used as a measure of accuracy and the required level of accuracy improves as the minimization proceeds and as (μa_k ) approaches a minimum. It is not necessary to compute the Hessian G as only the product Gd is required. The product is calculated using the finite difference formula given below:

where $\sigma =\frac{\sqrt{\mathit{machine\; precision}}}{\parallel \mathbf{d}\parallel }$. This avoids the calculation and storage of the Hessian but requires an additional gradient evaluation for each minor iteration. The pseudo codes for truncated Newton optimization is as follows:

The truncated Newton search direction is constructed so that it always satisfies Wolfe’s condition 1 (Step 3, ${\mathbf{d}}_k^T$ g _k ≤ -ε₁∥∥d _k ∥d _k∥). Wolfe’s conditions 2 and 3 are satisfied within the line search, as discussed is Step 4. Since all the Wolfe’s conditions have been satisfied, the method is globally convergent, and will terminate in a finite number of iterations.

A similar formula applies for lifetime τ and absorption coefficient μ_ax→m.

3.2 Reverse automatic differentiation for computation of ∇E_x,m

The reverse automatic differentiation (RAD) method is used to calculate the gradient g = g(μ̅_a or τ̅ or μ̅_ax→m). The reader interested in background literature is referred to References [28, 33, and 34]. A brief description of the reverse differentiation and analytic implementation of this method is discussed in the following section.

Reverse automatic differentiation is discussed in terms of the Wengert [35] list. This list decomposes complicated function of many variables into a sequence of simpler functions of one or two variables. Functions of two variables are called the “binary function”. Binary functions are either addition, or subtraction or multiplication. Functions of one variable are either reciprocation, raise of power, exponential, logarithm, trigonometric, or hyperbolic. These functions are defined as unary functions.

If f(x) is a function of the n variables x ₁,…,x_n then a set of new variables x _n+1,…x _P are introduced, where P is the number of arithmetic operations involved in calculating the function f(x). A Wengert list for the calculation for the calculation of f(x) consists of a list of these unary and binary function processed in a given order;

Reverse differentiation can be explained as follows. Given the Wengert list for the calculation of f a set of adjoint variables x̂ = ∂f /∂x_i is introduced and the order of the operations reversed. The method is based on the function of a function rule of calculus that implies

The automatic method involves a pass in the reverse direction through the list after the function has been evaluated in a forward sweep. Hence the method is:

Griewank [14] showed that it is always possible to calculate any gradient vector of a function in less than 3 times the computational cost of the function by the reverse automatic differentiation method irrespective of the dimension of the problem. The main disadvantage of the automatic code is that it requires large storage and significant data accessing overheads [28, 33, and 34]. In order to overcome these difficulties, the reverse differentiation method is extended and has been performed analytically. This does not require large storage and at the same time gradients are calculated in less than three times the operation count of the function. The gradients of E_x,m with respect to parameter μ_axf are calculated using the reverse differentiation method as follows.

The above gradients of error function E_x,m are calculated using the global stiffness matrices, K, and the matrix b assembled from the excitation and emission diffusion equations as described in Section 2.2 and 2.3. (The reverse differentiation method described above uses the fact that if Φ _x,m is the solution of the equation KΦ̅ _x,m = b then K̂ = ${-}_v^{\mathit{-}T}$ _x,m Φ̂ _x,m and b = v̅_x,m where v̅_x,m is the solution of the equation $\mathbf{K}{\overline{v}}_{x,m}={\hat{\overline{\Phi }}}_{x,m}$, see Appendix A). With ∇E_x,m and E_x,m computed, the Newton directions can be computed as described in Section 3.1. Similar formulas hold for τ̂ and μ̂_ax→m

4 Conclusions

We have formulated the solution to the absorption and fluorescence FDPM optical imaging problem as a non-linear optimization scheme. We develop the truncated Newton method for minimizing the sum of errors between measured and synthetic data generated using a simplistic finite element model with a minimal number of sources and detectors. Since the truncated Newton method requires only the product of the Hessian with a direction at any point, it is computationally efficient for large optimization problems. We couple the truncated Newton method with a finite element solver which assumes a polynomial series in absorption and lifetime in order to accurately compute Φ_x,m and employ principles of reverse differentiation in order to efficiently compute the gradient of the error function, ∇E_x,m . Using synthetic data with simulated measurement noise, we demonstrate in Part II the optimization strategy on a series of studies involving a minimum number of sources and detectors.