Fast tomographic reconstruction strategy for diffuse optical tomography

Yuhu Zhai; Steven A Cummer

doi:10.1364/OE.17.005285

1. Introduction

The technology development of Diffuse Optical Tomography (DOT) has been growing significantly in the past two decades as DOT is a promising tool for in vivo and non-invasive imaging of tissues using near-infrared light (NIR) [1–5]. DOT employs methodologies based on differences in the absorption and scattering properties between an “embedded object” such as cancer tumors and its surroundings. The differences in these optical properties between healthy and pathological tissues provide the contrast for diffuse optical imaging. Diffuse optical imaging to provide information at the molecular level is becoming one of the most promising alternative strategies to significantly advance biomedical research. This strategic approach could provide convenient, non-invasive imaging of tissues for early cancer detection and real-time tumor assessment [1–2]. It can also improve our ability to probe complex biologic interactions dynamically and to study disease and treatment responses over time.

The potential use of DOT and its capability to provide molecular imaging, however, are currently limited by excessively long image reconstruction time. Developing fast and robust reconstruction methods is the main challenge in making DOT a viable tool for clinical diagnostics [6]. Most current DOT software tools employ finite element (FE) models combined with appropriate tomographic principles to recover tissue optical properties from boundary measurements of light transmission [7–8], and to present them as 3-D images. For example, the TOAST software [9], developed by the Biomedical Optics Research Laboratory at the University College London, has been used as a time-resolved image reconstruction tool in experimental and clinical applications of DOT for neonatal imaging, breast imaging and muscle imaging. The NIRFAST software, developed by the Biolight research group at Dartmouth College [10], is a frequency-domain DOT imaging tool for breast cancer detection using near-infrared light. Even though such developed software can significantly reduce image reconstruction time compared to the MC model, it may still take hours to solve the imaging problem and present 3-D images to users [11–12]. This delay limits time-dependent applications such as the monitoring of cancer therapy, which would be possible only through a significant improvement in imaging speed. Recent clinical trials for diffuse optical techniques in breast imaging also show that imaging speed is the primary concern of patients [13].

It is important to develop faster solution strategies for real-time DOT imaging. Real-time imaging can provide important clinical benefits such as detection of cancer treatment needs at an early, relatively harmless stage [4]. The need for fast forward and inversion algorithms is becoming ever more important as measurement data sets increase in size due to larger and more precise NIR emitting arrays [1].

It has been demonstrated that the use of fast solvers [14] or the acceleration of numerical solutions using, for example, the adaptive mesh-based enhancement schemes [15] and the multi-grid methods [16] have been important contributions to achieving practical DOT schemes. However, the dual meshing scheme and multi-grid methods require carefully designed special treatment at mesh interfaces between healthy and pathological tissues depending on problems of interest. One potentially useful approach to overcome this computational difficulty is model order reduction (MOR), where the original large-order models obtained after FE discretization of the governing equations are replaced by models of substantially lower order, yet capable of capturing the behavior of the original subsystem with sufficient accuracy. By incorporating advanced MOR techniques, we introduce new imaging strategies that can demonstrate significantly reduced imaging time while preserving to the greatest extent tissue optical properties. The MOR-based imaging approach integrated with existing and new measurement systems for DOT will provide a variety of image reconstruction and data reduction techniques and close to real-time image reconstruction for disease screening. This is the first step toward development of imaging tools that will enable inexpensive, noninvasive, portable imaging systems that can run at patient’s bedside in real time and produce data continuously over a long period of time. In this research, we focus on the development of computationally efficient tomographic imaging strategies to explore fast imaging capability of DOT for clinical usage.

2. Imaging methods

Tissue optical imaging toward molecular level is the wave of the future. However, the long computational time of current model-based iterative algorithm may prevent DOT from providing molecular imaging capability. MOR is a numerical technique to extract essential characteristics of a large complex system in the form of a much smaller but sufficiently accurate model description. We develop computationally more efficient imaging tools based on the advanced MOR techniques to enable a fast DOT that could meet current clinical needs more effectively. Efficiency of the forward problem in the FE solver is critical for real-time imaging. The MOR-based imaging tool will be significantly faster (greater than 10 times) than Non-MOR based imaging tools, while being reasonably accurate. A “variable accuracy setting” to be controlled by users would be incorporated, for example, a user defined tolerance as an input to the software. The basic principle for the fast-solver technology lies in the advanced MOR techniques widely used in structural dynamics, control theory, VLSI (vary-large-scale-integrated) circuit simulations and MEMS (micro-electro-mechanical-system) simulations. The key idea underlying MOR-based image reconstruction is that instead of solving a large full-order model in the forward problem, which may have tens of thousands of unknowns, we will solve a reduced-order model, which is significantly smaller and faster to run than the large full-order model. Any existing fast FE-based numerical solvers, for example, the adaptive mesh-based enhancement schemes [15] and the multi-grid methods [16] as valuable contributions to achieving practical imaging schemes can be combined with the MOR techniques and thus may be implemented into the imaging tools as needed for specific applications. Fig. 1 shows the concept of the projection-based MOR. First, spatial discretization of the photon diffusion equation will result in a large full-order system to be solved. Three essential steps for MOR are then followed to reduce the large full-order system to a much smaller system that can be easily solved. The original system solution is obtained by back-projecting the unknowns from the small system to the large system.

Fig. 1. Concept of projection-based MOR. This figure depicts three essential steps to extract a reduced-order model from a large full-order model system through orthogonal projection.

Download Full Size | PDF

2.1. Forward problem and inverse problem

The finite element method has been proven to be powerful in studies of both forward and inverse photon migration problems. Based on the diffusion equation that describes the propagation of light in multiply scattering media, the FE-based reconstruction algorithm is centered on Newton's iterative method for minimizing the error between FE solution and measurement data on the object boundaries [8]. For iteration, the FE solution is followed by a synthesized scheme of Marquardt and Tikhonov regularizations to obtain the updated tissue parameter distribution [2, 17–18].

The diffusion approximations to the radiative transport equation [19] in the time-domain and frequency-domain are

\frac{1}{c} \frac{\partial ρ (r, t)}{\partial t} + (μ_{a} (r)) - \nabla \cdot (κ (r) \nabla)) ρ (r, t) = s (r, t),

where ρ(r,t) or ρ(r,ω) is the photon density, s the source strength, μ_a (r)the absorption coefficient, κ(r) the diffusion coefficient, ω the excitation frequency and c the speed of light in the medium. The FE-based approaches recast the image reconstruction problem into a nonlinear optimization problem, where the optimization parameters are coefficients in a basis function expansion for the spatial distribution of tissue optical properties. Specifically, let $ρ (r, t) = \sum_{i = 1}^{n} N_{i} (r) u_{i} (t)$ be a n-node FE discretization of the photon density, N_i the interpolation function, r the position vector in a domain of imaging Ω with ∂Ω the domain boundary, and u = [u ₁(t) u ₂(t) ⋯u_n(t)]^T the nodal photon density. The result after Galerkin FE discretization can be expressed in matrix form

M \frac{d u}{d t} + Ku = s,

where M, K = B(κ) + C(μ_a), s are system matrices with

B_{ij} = \int_{Ω} κ (r) \nabla N_{j} (r) \cdot \nabla N_{i} (r) d Ω; C_{ij} = \int_{Ω} μ_{a} (r) N_{j} (r) N_{i} (r) d Ω;

M_{ij} = \int_{Ω} N_{j} (r) N_{i} (r) dΩ; s_{j} (t) = \int_{Ω} N_{j} (r) {S (r, t) d Ω + \int_{\partial Ω} N}_{j} (r) Γ (r, t) d \partial Ω;

in time domain, and

Ku = s, C_{ij} = \int_{Ω} (μ_{a} (r) + \frac{iω}{c}) N_{j} (r) N_{i} (r) d Ω

in frequency domain. We then solve the inverse problem by performing iteration of the trial images based on assumed tissue optical properties. The solution is obtained by minimization of the difference error of the data from output of the forward model on such trial solutions.

The inverse problem involves minimizing the objective function defined as

χ^{2} = ∥ {\frac{\partial u}{\partial θ} ∥ Δ θ - Δ θ}^{2} + λ {∥ Δ θ ∥}^{2}

to simultaneously recover optical properties μ_a(r) and κ(r) in group form θ = (μ_a,κ). Minimization of the objective function provides the update equation for iterations

Δ θ = {(J^{T} J + λI)}^{- 1} J^{T} Δ u,

where λ is the regularization parameter and $J = \frac{\partial u}{\partial θ}$ is the Jacobian matrix.

Although the FE-based iterative approach is powerful, it requires a considerable amount of time and computational power due to nonlinearity of tissue optical properties [11]. For example, for a forward problem solution, a mesh size of over 1000 elements may be necessary [2] and the nonlinear iterations for 2-D imaging may take an hour to run on a Pentium PC [15]. Therefore, for DOT to realize its potential for real-time imaging, an efficient software tool is needed in refining the forward solution scheme for its appropriate integration with an experimental data acquisition system.

2.2. Model-order reduction

We introduce a MOR-based forward modeling strategy and a robust de-convolution technique for fast tomographic inversion. The MOR-based modeling technique provides means to significantly reduce imaging time while maintaining reasonable accuracy. The basic idea is to find a convenient and inexpensive transformation of coordinates to simplify the governing equations that are easier to solve. Conventional mode superposition and the component-mode synthesis methods for model-order reduction require solving a large and expensive eigenvalue problem [20]. The use of eigenvectors, to reduce the size of systems or to represent the system behavior by a small number of generalized coordinates, is not necessarily the best approach. Several alternative methods for model reduction without solving an eigenvalue problem are the Lanczos method and the WYD (Wilson-Yuan-Dickens) method, belonging to the well-known Krylov-subspace methods [20–21]. The Krylov-type MOR methods generate a set of basis vectors by repetitively multiplying a starting vector by the spatially discretized system matrices. The system information is well captured by a subspace spanned by the Krylov basis vectors generated at a small fraction of overall computational cost.

For a nonlinear iteration, instead of solving a full-order forward model, we solve a reduced-order model with a dimension r by invoking the following coordinate transformation

u_{n \times 1} = Q_{n \times r} η_{r \times 1,}

where η = [η ₁(t) η ₂(t)⋯η_r(t)]^T is the discrete photon density in the transformed coordinate. The transformation matrix Q consists of WYD basis vectors in the Krylov-subspace κ ^(r)(K ^-1 M,w ₁) = span[w ₁ (K ^-1 M)w ₁ (K ^-1 M)² w ₁ ⋯ (K ^-1 M)^r-1 w ₁] with the starting vector w ₁ = K ^-1 s the static solution of the problem; or the Lanczos basis vectors in the Krylov-subspace κ ^(r)(K ^-1 M,z ₁) = span[z ₁ (K ^-1 M)z ₁ (K ^-1 M)² z ₁ ⋯ (K ^-1 M)^r-1 z ₁] where a random starting vector is used. The WYD algorithm is implemented in the image reconstruction to generate trial vectors. The algorithms for generating Lanczos basis vectors and the WYD basis vectors are presented. The generalized Falk method [20] can be viewed as an inverse WYD method, where the Krylov sequence is generated from an inverse matrix of (K ^-1 M). The Pade-via-Lanczos (PVL) method, widely used in VLSI circuit simulations, is the frequency-domain Lanczos method based on Pade approximation through an effective method of moment-matching [21].

2.2.1. Lanczos method and WYD method

The Lanczos method was originally proposed as a technique for matrix tridiagonalization. A sequence of trial vectors are formed by repeatedly multiplying the matrix to be reduced with a random starting vector and these vectors form a Krylov sequence. Each new vector is orthogonalized with respect to the two previous vectors. This orthogonalization procedure can be shown to be sufficient to obtain orthogonality with all previously calculated vectors. The coefficients computed from the orthogonalization process are then combined to form a tridiagonal matrix that theoretically has the same eigenvalues as the original matrix, after all n Lanczos vectors are generated.

Lanczos algorithm

Given M, K ∈ R ^n×n and s ∈ R ^n×1
Triangularize matrix K = LDL ^T
Obtain starting vector z ₁ by normalizing a randomly generated vector z ^* ₁ with respect to M such that $z_{1} = \frac{z_{1}^{*}}{\sqrt{z_{1}^{*^{T}} M z_{1}^{*}}}$
Solve for additional vectors with b ₁ = 0 and i = 2,…, r, with r << n $K z_{i}^{*} = M z_{i - 1} \Rightarrow z_{i}^{*} = K^{- 1} M z_{i - 1}$
Construct Lanczos Ritz vector matrix for coordination transformation Z = [z ₁ z ₂ … z _r]_n×r
Construct symmetric tridiagonal matrix of reduced order r, $T_{r} = {[\begin{matrix} c_{1} & b_{2} & 0 & \dots & \dots & 0 \\ b_{2} & c_{2} & b_{3} & \dots & \dots & 0 \\ 0 & b_{3} & c_{3} & \dots & \dots & 0 \\ 0 & \dots & \dots & \dots & \dots & 0 \\ 0 & \dots & \dots & b_{r - 1} & c_{r - 1} & b_{r} \\ 0 & \dots & \dots & 0 & b_{r} & c_{r} \end{matrix}]}_{r \times r} .$

Although similar to Lanczos method in generating the Krylov sequence, the WYD method is based on the direct superposition of Ritz vectors constructed from the spatial distribution of the specific dynamic excitation, which is important information neglected by direct use of exact mode shapes in conventional mode superposition. Unlike the Lanczos method, the first WYD vector is the static response obtained from the spatial distribution of the dynamic excitation. It ensures that the generated trial vectors contain important information of the system dynamic response.

WYD algorithm

Given M, K ∈ R ^n×n and s ∈ R ^n×1
Triangularize matrix K = LDL ^T
Solve for excitation-dependent starting vector w ^* ₁ = K ^-1 s, and normalize w ^* ₁ to obtain $w_{1} = \frac{w_{1}^{*}}{\sqrt{w_{1}^{*^{T}} M w_{1}^{*}}}$
Solve for additional vectors i = 2,…,r,
$K w_{i}^{*} = M w_{i - 1} \Rightarrow w_{i}^{*} = K^{- 1} M w_{i - 1}$
Construct WYD Ritz vector matrix for coordination transformation W = [w ₁ w ₂ … w _r]_n×r

The basic algorithm to generate the WYD vectors is equivalent to the Lanczos method. Both methods can be used to reduce the size of the original model by truncating the generation of the Ritz vectors to obtain an approximation to the original eigenpairs. In the WYD method, each new Ritz vector is generated by an orthogonalization with all previously generated vectors; in the Lanczos method, each new vector is generated by an orthogonalization with only the two previously generated vectors. If we write the orthogonalization step in WYD algorithm in matrix form, we obtain

K^{- 1} M [w_{1} w_{2} \dots w_{r}] = [w_{1} w_{2} \dots w_{r}] H,

where entries of matrix

H = {[\begin{matrix} c_{2,1} & c_{3,1} & c_{4,1} & \dots & \dots & c_{r + 1,1} \\ b_{2} & c_{3,2} & c_{4,2} & \dots & \dots & c_{r + 1,2} \\ 0 & b_{3} & c_{4,3} & \dots & \dots & \dots \\ 0 & 0 & b_{4} & \dots & \dots & \dots \\ 0 & \dots & \dots & b_{r - 1} & c_{r, r - 1} & c_{r + 1, r - 1} \\ 0 & \dots & \dots & 0 & b_{r} & c_{r + 1, r} \end{matrix}]}_{r \times r}

in the upper Hessenberg form consist of orthogonalization coefficients. Premultiplying Eq. (9) by W ^T M, and use the fact that the technique used to construct the WYD basis vectors enforces the orthonormality, i.e., W ^T MW = I ∈ R ^r×r, we obtain

W^{T} M K^{- 1} MW = H,

The left hand side of Eq. (10) must be symmetric since K and M are symmetric. The only way that Eq. (10) can be satisfied is when H is tridiagonal such that H = T _r. Therefore, the generation of WYD basis vectors is equivalent to the generation of the Lanczos vectors.

Substituting Lanczos coordinate-transformation u = Zη into the discrete system equations in Eq. (2) and premultiplying the result first by MK ^-1, then by Z ^T, we obtain

(Z^{T} M K^{- 1} MZ) \frac{d η}{dt} + (Z^{T} MZ) η = Z^{T} M K^{- 1} s .

Using the fact that H = T _r and Z ^T MZ = I ∈ R ^r×r, we obtain from Eqs. (10) and (11)

T_{r} \frac{d η}{dt} + Iη = Z^{T} M K^{- 1} s .

Similarly, substituting WYD coordinate-transformation u = Wη into the discrete system equations in Eq. (2) and premultiplying the result by W ^T, we obtain

I \frac{d η}{dt} + K^{*} η = W^{T} s .

The reduced matrix K = W ^T KW ∈ R ^r×r is in general a full matrix. The reduced-order model in Eq. (12) or Eq. (13) is a set of r equations that can be solved by a direct step-by-step integration method, or by introducing an additional transformation for a reduction to a diagonal form (mode superposition method). We use the explicit time-integration scheme (such as the Newmark algorithm) to solve the reduced system equations to avoid solving the reduced order eigenvalue problem. Once solution of the reduced order system is obtained, solution in the original physical coordinate is easily recovered by matrix multiplication as shown in Eq. (8), with Q = Z for Lanczos method or Q = W for WYD method.

A frequency-domain reduced-order model is obtained by projecting the original system Eq. (5) onto the subspace spanned by the WYD basis vectors as below

K^{*} η = s^{*} (η),

where K ^* = B ^* + C ^*, and B ^* = W ^T BW ; C ^* = W ^T CW; s ^* = W ^T s are the system matrices in the reduced coordinates with much smaller dimensions. Accuracy of the reduced-order model can be controlled by a user defined tolerance e based on the truncation of projected load s ^*as compared with the original input. For a given e, this truncation based on e = ∥s - Us ^*∥ will automatically generate the order of the reduced-order model.

We implement the above Krylov-type methods into the frequency-domain DOT imaging procedure. The overhead computation time for generating a few trial vectors in the Krylov sequence is negligible [20–22]. A frequency-domain model is useful for reconstruction of tissue internal optical property distributions for each wavelength. A time-domain model is necessary when phosphorescence or inelastic scattering within the tissue cannot be neglected and a change in tissue optical properties can be captured in time. The MOR-based forward solution procedure using Krylov-type methods is particularly useful for weakly non-linear optical imaging to find the difference of optical images where measurements are taken before and after a small change in tissue optical properties [23]. In such cases, the actual optical properties are close to an initial estimate and the measured data are close to the simulated measurements. Fig. 2 shows the procedure of an MOR-based image reconstruction for DOT. Although the geometry of tissue structures does not permit a discretization within a few finite elements, the photon behavior in tissues may be well characterized by a few generalized degrees of freedom.

Fig. 2. The MOR-based tomographic reconstruction.

Download Full Size | PDF

2.2.2. Proper orthogonal decomposition

The time-resolved image reconstruction will employ nonlinear MOR techniques such as the proper orthogonal decomposition (POD) [22], or the trajectory piecewise-linear approach [24] for efficient DOT imaging. We implement the fixed-point solution technique, combined with the POD-based MOR technique to efficiently solve the forward problem. The idea of fixed-point solution technique is to shift the nonlinear dependence of tissue optical properties to the right hand side of the system equation. The fixed-point technique is particularly suitable for POD based model-order reduction because the basis vectors in the transformation matrix do not need to be updated for nonlinear iterations. Only the right hand side term of the equation needs to be updated. Newton's method involves the derivation of a tangent stiffness matrix by a linearization of the nonlinear term and thus requires a more careful formulation. Moreover, the fixed-point technique is easier to implement and it has been successfully used for POD-based MOR to solve computational electromagnetic problems with nonlinear hysteresis in power electronics [25–26]. The POD-based reconstruction technique has also been used successfully for radio tomographic imaging of Earth's magnetosphere [27].

The idea of POD is to extract a typical spatial structure from an ensemble of realization (called snapshots) of a dynamic system. In the proposed method, the basis functions for coordinate transformation may be generated numerically from a small number of “snapshots”, of the state variables sampled in time. In a finite element discretization of space variable such as photon density at n discrete spatial points (nodes), the sample distributions U = [u ₁ u ₂ ⋯ u _r]_n×r are a series of snapshots or time frames. The empirical basis vectors can be selected as the normalized left singular vectors in the singular value decomposition (SVD) of the rectangular snapshot matrix U,

U = Ψ \sum V^{T} .

We use the reduced SVD of U to efficiently compute the POD basis vectors. Once constructed, the reduced order model can be executed repeatedly for problems with a variety of different excitations at very low computational cost. It is also important to take enough snapshots when the dynamic character of the system varies rapidly.

Instead of solving the full-order model system Eq. (2), which may involve thousands of state variables, we express the unknowns in the subspace spanned by the POD basis vectors

u_{n \times 1} = Ψ_{n \times r} η_{r \times 1},

where r << n. We project the semi-discrete Eq. (2) onto the POD basis to obtain the reduced-order system equation with much smaller dimensions, such that

M^{*} \frac{d η}{dt} + K^{*} η = s^{*},

where M ^* = Ψ ^T MΨ, K ^* = B ^* + C ^*, and B ^* = Ψ ^T BΨ; C ^* = Ψ ^T CΨ ; s ^* = Ψ ^T s are system matrices in the subspace spanned by the POD basis vectors. The explicit time-integration scheme is used to obtain solution of the reduced-order system and solution in the physical coordinate is easily recovered by matrix multiplication as shown in Eq. (16 ).

The software tool can be integrated into the DOT medical imaging system. It is also straightforward to incorporate statistical-based a priori information into the MOR-based imaging scheme, as we have implemented successfully elsewhere [27].

3. Numerical results

We have developed and successfully implemented linear and nonlinear MOR techniques into the time-domain and frequency-domain finite element forward solver for the selection of the best MOR strategy for DOT imaging. Fig. 3 presents the simulated photon density distribution of a test phantom [8] from the full-order model and that from the reduced-order model. The reduced-order model is generated using the WYD method described earlier. The reduced-order model is over 10 times faster while maintaining a 5% overall mean square error norm. Fig. 4 presents distribution of the original absorption and reduced scattering coefficients in the test phantom. The background absorption coefficient is μ_a = 25 m^-1, the reduced scattering coefficient is μ′_s = 2000 m^-1, the diffusion coefficient $κ = \frac{1}{3 (μ_{a} + μ_{s}^{'})} = 1 / 6075 m$ ; the target absorption coefficient is μ_a = 50 m^-1 and the target diffusion coefficient κ = 1/12150 m. Fig. 5 presents the reconstructed absorption coefficient and reduced scattering coefficient.

Fig. 3. The simulated photon density distribution of a test phantom from the full-order model (left) and that from the reduced-order model (right).

Download Full Size | PDF

Fig. 4. The original absorption μ_a (left) and reduced scattering μ′_s (right) coefficients of a test phantom.

Download Full Size | PDF

Fig. 5. The reconstructed absorption μ_a (left) and reduced scattering μ′_s (right) coefficients of a test phantom.

Download Full Size | PDF

Efficiency of the proposed imaging tools depends on the MOR techniques used and their feasibility and compatibility with existing FE-based DOT imaging procedures. We found that the WYD method is more advantageous in terms of efficiency and controllability of reconstruction error as compared to the generalized Falk method and the PVL method for a weakly nonlinear DOT image reconstruction (see Fig. 6). However, the results from our preliminary study demonstrate that the proposed MOR-strategy has the potential to become an effective means to significantly reduce imaging time to achieve real-time optical imaging.

Table 1 shows the comparison of computational speed for a time-dependent FE simulation using the full-order model and the reduced-order model for the forward problem of DOT. The FE mesh has a dimension of 282 triangular elements and 159 nodes. The resulting full-order model with 159 unknowns is reduced to 40 unknowns in the reduced-order model using the WYD method. The speed-up ratio from the MATLAB solver is 12 for 200 time steps. The overhead cost for generating the basis vectors is 0.17s, which is only 1/10 of the full-order simulation time. For a refined mesh with 599 nodes, the speed-up ratio is 33. The quality of the reconstructed image depends strongly on accuracy of the forward model. The mean square error between the full-order model solution and the reduced-order model solution for photon density at the final time step is less than 1%. The left panel of Fig. 6 shows the error norm of the simulated photon density between the full-order model and the reduced-order model. It clearly demonstrated that the reduced-order model converges to the full-order model with increased number of WYD trial vectors. The reduced-order model, however, is significantly fast to run than the full-order model. The right panel of Fig. 6 shows the error norms for various Krylov-type MOR methods for a forward diffusion problem. The WYD method is clearly more advantageous than the generalized Falk method and the PVL method in terms of convergence and it will thus be selected for the MOR-based imaging software development.

Table 1. Speed-up ratio for linear forward problem of DOT using the WYD method

View Table | View all tables in this article

Fig. 6. The reduced-order model converges to the full-order model with increased number of trial vectors for MOR and comparison of various Krylov-type MOR methods.

Download Full Size | PDF

Table 2 shows the speed comparison for a time-dependent FE simulation between the full-order model and the reduced-order model for the forward and the inverse problem in DOT imaging. The speed-up ratio from a MATLAB software tool is 3.5 for 20 time steps. Due to long simulation time for the forward and inverse combined problem in DOT imaging, we only run 20 time steps for speed comparison. However, this result is a very important proof that the WYD method works reasonably well for the nonlinear DOT imaging.

Table 2. Speed-up ratio for nonlinear DOT imaging using the WYD method

View Table | View all tables in this article

We have also implemented the POD method for efficient DOT imaging based on the Newton iterative reconstruction procedure using the same FE mesh. For a time-dependent simulation with 40 steps, the POD-based imaging reconstruction is 11 times faster than the reconstruction without MOR. The basis vectors are generated from the SVD of 40 snapshots of a full-order simulation with similar but different source excitations. The error tolerance defined by the mean square error between 42 boundary measurements and the FE calculated photon density is 1×10^-8 for the full-order model, and 5×10^-8 for the reduced-order model.

Table 3. Speed-up ratio for nonlinear DOT imaging using the POD-based MOR technique

View Table | View all tables in this article

The MOR-based DOT techniques have demonstrated reasonably accurate but significantly fast DOT imaging capabilities. The numerical results provide a very important proof of the proposed MOR concept for DOT imaging. The developed image reconstruction tool can be installed in a medical imaging system for DOT, such as the Laser Breast Scanner or a portable optical scanner for early detection of breast cancer. The software tool has high potential to be used for real-time optical imaging and image data analysis in conjunction with experimental systems such as the NIR-Scanner. The fast imaging tools are developed based on Krylov-type MOR techniques such as the WYD method, and the nonlinear MOR techniques such as the POD method implemented in our study. The frequency-domain solver will be used for spectral analysis and the time-domain solver will be used for applications that require real-time monitoring of response to therapy. The main goal is to build up a platform that can be used to develop a real-time imaging and image data analysis system.

6. Summary

We have developed fast imaging strategies for solving the forward problem and the inverse problem by implementing two Krylov-type MOR techniques into the FE-based imaging procedure using MATLAB. This preliminary software prototype has been used to demonstrate efficiency of the proposed imaging technique that can be integrated with an experimental data acquisition system. Our numerical study shows that the WYD method is more efficient compared to the generalized Falk method, and more accurate compared to the PVL method. This is mainly because lower-order modes in an optical diffusion system are better captured by the WYD trial vectors, as compared to the generalized Falk method, where higher order modes are better described. A selective reorthogonalization technique [22] is used to generate the WYD trial vectors, while no reorthogonalization is used in the PVL method. The software tool developed and tested can be translated into FORTRAN or C for higher efficiency and greater flexibility for its integration with experimental data acquisition system. The ultimate goal is to achieve real-time imaging. This can be demonstrated with true experimental data in future work. The experimental data can be obtained through the breast cancer multidimensional diffuse optical imaging network [28–29].

Due to its linear nature for characterizing a nonlinear dynamic system, the POD-based MOR technique has become popular in the past few years and has been widely used for nonlinear MOR in various engineering applications. We have demonstrated in our preliminary study that the POD method can greatly improve the imaging reconstruction speed for DOT. The POD technique for MOR has the great potential for time-dependent DOT imaging problem, which currently may take an hour for 3-D imaging [11], and efficiency is currently the bottleneck that may prevent the clinical use of DOT [6]. We have successfully applied the POD technique to solve radio tomographic imaging problems when combined with the robust de-convolution technique [27]. We have also incorporated POD technique into the MOR in the preliminary study, which showed much faster computational speed than full-order model simulation.

7. Conclusion

We have developed a MOR-based image reconstruction technique to significantly improve computational efficiency of DOT. We have implemented advanced linear and nonlinear MOR techniques for reducing system complexity, while preserving (to the greatest possible extent) system input-output behavior for the forward problem. The nonlinear DOT imaging using WYD MOR method is significantly faster than the DOT imaging without MOR. The nonlinear DOT using POD-based MOR method is an order of magnitude faster than that without using the POD method.

References and links

1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nature Biotechnol. 23, 313–320 (2005). [CrossRef]

2. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R01–R43 (2005). [CrossRef]

3. B.W. Pogue, S.C. Davis, X. Song, B. A. Brooksby, H. Dehghani, and K. D. Paulsen, “Image analysis methods for diffuse optical tomography,” J. Biomed. Opt. 11, 033001 (2006). [CrossRef]

4. S.L. Jacques and B.W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt. 13, 041302 (2008). [CrossRef] [PubMed]

5. B. J. Tromberg, B. W. Pogue, K.D. Paulsen, A.G. Yodh, D.A. Boas, and A.E. Cerussi, “Assessing the future of diffuse optical imaging technologies for breast cancer management,” Med. Phys. 35, 2443–2451 (2008). [CrossRef] [PubMed]

6. M. Schweiger, A. Gibson, and S. R. Arridge, “Computational aspects of diffuse optical tomography,” Comput. Sci. Eng. 5, 33–41 (2003). [CrossRef]

7. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

8. K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys. 22, 691–701 (1995). [CrossRef] [PubMed]

9. The Biomedical Optics Research Laboratory at University College London, http://www.medphys.ucl.ac.uk/research/borl/.

10. The Near Infrared Imaging Group at Dartmouth College, http://www-nml.dartmouth.edu/nir/.

11. P. K. Yalavarthy, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Critically computational aspects of near infrared circular tomographic imaging: Analysis of measurement number, mesh resolution and reconstruction basis,” Opt. Express 14, 6113, (2006). [CrossRef] [PubMed]

12. P. K. Yalavarthy, D.R. Lynch, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Implementation of a computationally efficient least-squares algorithm for highly under-determined three-dimensional diffuse optical tomography problems,” Med. Phys. 35(5), 1682–1696, (2008). [CrossRef] [PubMed]

13. N. S. Shah, A. E. Cerussi, D. Gordon, A. Durkin, L. Wenzel, B. Hill, M. Compton, and B. J. Tromberg, “Integration of diffuse optical technology into clinical settings for breast health applications,” Frontiers in Optics, The 90^th OSA Annual Meeting, 2006.

14. J. Ripoll, M. Nieto-Vesperinas, R. Weissleder, and V. Ntziachristos, “Fast analytical approximation for arbitrary geometries in diffuse optical tomography,” Opt. Lett. 27, 527, (2002). [CrossRef]

15. X. Gu, Y. Xu, and H. Jiang, “Mesh-based enhancement schemes in diffuse optical tomography,” Med. Phys. 30, 861–869 (2003). [CrossRef] [PubMed]

16. J. C. Ye, C. A. Bouman, K. J. Webb, and R. P. Millane, “Nonlinear multigrid algorithms for Bayesian optical diffusion tomography,” IEEE Trans. Image Process. 10, 909–922 (2001). [CrossRef]

17. H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom and clinical results,” Appl. Opt. 42, 135–145 (2003). [CrossRef] [PubMed]

18. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys. 34(6), 2085–2098 (2007). [CrossRef] [PubMed]

19. S. R. Arridge and J. C. Hebden. “Optical imaging in medicine: II. Modeling and reconstruction,” Phys. in Med. Biol. 42, 841–853 (1997). [CrossRef]

20. L. Vu-Quoc, Y. Zhai, and K. D. T. Ngo, “Efficient simulation of coupled circuit-field problems: Generalized Falk method,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. 23, 1209–1219 (2004). [CrossRef]

21. P. Feldmann and R. Freund, “Reduced order modeling of large linear subcircuits via a block Lanczos algorithm,” In Proceedings Design Automation Conference (Piscataway, NJ, USA, 1995), pp.474–479.

22. Y. Zhai, “Model-order reduction for efficient simulation of nonlinear electro-magneto-thermal coupled problems,” Ph.D. Thesis, University of Florida, 2003.

23. J. Phillips, “Projection-based approaches for model reduction of weakly nonlinear, time-varying systems,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. 22, 171–187 (2003). [CrossRef]

24. M. Rewienski and J. White, “A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. 22, 155–169 (2003). [CrossRef]

25. Y. Zhai and L. Vu-Quoc, “Analysis of power magnetic components with nonlinear static hysteresis: Finite-element formulation,” IEEE Transactions on Magnetics 41, 2243–2256 (2005). [CrossRef]

26. Y. Zhai and L. Vu-Quoc, “Analysis of power magnetic components with nonlinear static hysteresis: Proper orthogonal decomposition and model reduction,” IEEE Trans. Magnetics 43, 1888–1897 (2007). [CrossRef]

27. Y. Zhai and S. A. Cummer, “An orthogonal projection and regularization technique for magnetospheric radio tomography,” J. Geophys. Res. 111, A03207 (2006). [CrossRef]

28. B. Brooksby, S. Jiang, H. Dehghani, B. W. Powgue, and K. D. Paulsen, “Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. of Biomed. Opt. 10, 0515041–10 (2005).

29. The Breast Cancer Multi-Dimentional Diffuse Optical Imaging Network, http://www.bli.uci.edu/ntroi/.

Full order	Reduced order	Speed-up ratio	Error norm
n=159	r=40	12.0	0.53%
n=599	r=40	32.6	0.89%

Full order	Reduced order	Speed-up ratio	Error norm
n=159	r=40	12.0	0.53%
n=599	r=40	32.6	0.89%

Fast tomographic reconstruction strategy for diffuse optical tomography

Abstract

1. Introduction

2. Imaging methods

2.1. Forward problem and inverse problem

2.2. Model-order reduction

2.2.1. Lanczos method and WYD method

2.2.2. Proper orthogonal decomposition

3. Numerical results

6. Summary

7. Conclusion

References and links

Cited By

Figures (6)

Tables (3)

Equations (30)

Optics Express