## High-performance image reconstruction in fluorescence tomography on desktop computers and graphics hardware |

Biomedical Optics Express, Vol. 2, Issue 11, pp. 3207-3222 (2011)

http://dx.doi.org/10.1364/BOE.2.003207

Acrobat PDF (971 KB)

### Abstract

Image reconstruction in fluorescence optical tomography is a three-dimensional nonlinear ill-posed problem governed by a system of partial differential equations. In this paper we demonstrate that a combination of state of the art numerical algorithms and a careful hardware optimized implementation allows to solve this large-scale inverse problem in a few seconds on standard desktop PCs with modern graphics hardware. In particular, we present methods to solve not only the forward but also the non-linear inverse problem by massively parallel programming on graphics processors. A comparison of optimized CPU and GPU implementations shows that the reconstruction can be accelerated by factors of about 15 through the use of the graphics hardware without compromising the accuracy in the reconstructed images.

© 2011 OSA

## 1. Introduction

1. S. Mordon, V. Maunoury, J. M. Devoisselle, Y. Abbas, and D. Coustaud, “Characterization of tumorous and normal tissue using a pH-sensitive fluorescence indicator (5,6-carboxyfluorescein) in vivo,” J. Photochem. Photobiol. B **13**, 307–314 (1992). [CrossRef] [PubMed]

2. I. Gannot, I. Ron, F. Hekmat, V. Chernomordik, and A. Gandjbakhche, “Functional optical detection based on pH dependent fluorescence lifetime,” Lasers Surg. Med. **35**, 342–348 (2004). [CrossRef] [PubMed]

3. E. Shives, Y. Xu, and H. Jiang, “Fluorescence lifetime tomography of turbid media based on an oxygen-sensitive dye,” Opt. Express **10**, 1557–1562 (2002). [PubMed]

4. B. Zhang, X. Yang, F. Yang, X. Yang, C. Qin, D. Han, X. Ma, K. Liu, and J. Tian, “The CUBLAS and CULA based GPU acceleration of adaptive finite element framework for bioluminescence tomography,” Opt. Express **18**, 20201–20214 (2010). [CrossRef] [PubMed]

5. J. Prakash, V. Chandrasekharan, V. Upendra, and P. K. Yalavarthy, “Accelerating frequency-domain diffuse optical tomographic image reconstruction using graphics processing units,” J. Biomed. Opt. **15**, 066009 (2010). [CrossRef]

*AGILE*[6] which is available at http://www.imt.tugraz.at/research/agile-gpu-image-reconstruction-library/.

## 2. Methods

### 2.1. Mathematical model

7. S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl. **15**, R41–R93 (1999). [CrossRef]

8. A. Joshi, W. Bangerth, and W. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express **12**, 5402–5417 (2004). [CrossRef] [PubMed]

*φ*,

_{i}*i*=

*x*,

*m*are the photon densities for the excitation (x) and the emission (m) light in the domain Ω covered by tissue, and the source term

*q*models a collimated light source. The parameters

*κ*and

_{i}*μ*are the diffusion and absorption coefficients of the tissue, and

_{i}*γ*is a measure incorporating the fluorophore’s extinction coefficient and its quantum yield. We assume that these parameters depend in a known form on the concentration

*c*of fluorophore in the tissue, and we write

*κ*(

_{i}*c*),

*μ*(

_{i}*c*) or

*γ*(

*c*) if this dependence needs to be emphasized.

7. S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl. **15**, R41–R93 (1999). [CrossRef]

*q*.

*d*allows to model the aperture and additional (linear) transfer characteristics of the detector. Other types of measurement functions are possible as well and have been used in e.g. [9

9. S. R. Arridge and M. Schweiger, “The use of multiple data types in time-resolved optical absorption and scattering tomography (TOAST),” in “Mathematical Methods in Medical Imaging II,”, D. C. Wilson and J. N. Wilson, eds. (1993), *Proc. SPIE*2035, pp. 218–229. [CrossRef]

11. S. R. Arridge, “Photon-measurement density functions. part I: Analytical forms,” Appl. Opt. **34**, 7395–7409 (1995). [CrossRef] [PubMed]

*q*,

_{j}*j*= 1,...,

*ns*, and taking measurements of the emission light with an array of detectors

*d*,

_{i}*i*= 1,...,

*nd*at the boundary. This gives rise to a measurement matrix ℳ of dimension

*nd*×

*ns*, which obviously depends on the distribution

*c*of the fluorophore, and we write ℳ(

*c*) to express this dependence.

### 2.2. Image reconstruction

*c*which gave rise to the measurements ℳ

*—which are perturbed by measurement noise and model errors—of the emitted light. A Tikhonov-type regularization term is incorporated in the objective functional for stability reasons. Then the reconstruction problem reads To solve this problem in a fast and stable manner, we utilize the iteratively regularized Gauß-Newton method (IRGN) [12], which leads to the iterative algorithm Here,*

^{δ}*c*onto the output ℳ(

*c*) and

*α*stabilizes the solution and is usually decreased in every iteration.

_{k}13. M. Hanke, “Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems,” Numer. Func. Anal. Optim. **18**, 971–993 (1997). [CrossRef]

14. D. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. **11**, 431–441 (1963). [CrossRef]

15. B. Kaltenbacher, “Some Newton-type methods for the regularization of nonlinear ill-posed problems,” Inv. Probl. **13**, 729–753 (1997). [CrossRef]

8. A. Joshi, W. Bangerth, and W. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express **12**, 5402–5417 (2004). [CrossRef] [PubMed]

16. H. Jiang, “Frequency-domain fluorescent diffusion tomography: A finite-element-based algorithm and simulations,” Appl. Opt. **37**, 5337–5343 (1998). [CrossRef]

19. H. Egger and M. Schlottbom, “Efficient reliable image reconstruction schemes for diffuse optical tomography,” Inv. Probl. Sci. Eng. **19**, 155–180 (2011). [CrossRef]

20. A. Godavarty, E. M. Sevick-Muraca, and M. J. Eppstein, “Three-dimensional fluorescence lifetime tomography,” Med. Phys. **32**, 992–1000 (2005). [CrossRef] [PubMed]

21. M. Freiberger, H. Egger, and H. Scharfetter, “Nonlinear inversion schemes for fluorescence optical tomography,” IEEE Trans. Biomed. Eng. **57**, 2723–2729 (2010). [CrossRef]

## 3. Implementation

### 3.1. Discretization

*N*vertices consisting of

_{V}*N*tetrahedral elements. For computational efficiency but also simplicity, the optical parameters and the fluorophore concentration

_{T}*c*are discretized by piecewise constant functions over the same mesh. The discretized version of the system (1)–(2) reads In this linear system,

*K*,

_{i}*M*and

_{i}*R*,

_{i}*i*∈ {

*x*,

*m*}, denote the stiffness-, mass- and boundary mass-matrices at the respective wavelength. Each column of the

*N*×

_{V}*ns*matrix

*Q*corresponds to a different excitation

*q*,

_{j}*j*= 1,...,

*ns*, and the

*ns*columns of the solution matrices

*V*and

_{x}*V*store the vector representations of the corresponding photon density fields

_{m}*φ*and

_{x}*φ*. The measurement data are given by ℳ(

_{m}*c*) =

*D*

^{⊤}

*V*, where each of the

_{m}*nd*columns of the detector matrix

*D*corresponds to one of the detectors

*d*,

_{i}*i*= 1,...,

*nd*, in (4).

*S*=

*S*(

*c*) is an

*nd*×

*ns*×

*N*tensor (respectively an (

_{T}*nd*·

*ns*) ×

*N*matrix). Its action onto a parameter perturbation

_{T}*h*∈ ℝ

*is given by where*

^{NT}*h*is an arbitrary element from the

*N*dimensional parameter space, and

_{T}*W*,

_{m}*W*denote the solutions of the adjoint system

_{x}### 3.2. Initialization

*l*and

*l*+ 1 are stored as sparse prolongation matrices

*l*to the next finer level

*l*+ 1, i.e.

*u*to the finer mesh where it is represented by

_{l}*u*

_{l}_{+1}. The mesh hierarchy is only required for the initialization step. The memory requirement for storing the prolongation matrices

*O*(

*N*).

_{V}*N*× 4 table containing the global vertex numbers for the

_{T}*N*individual elements. The 4-element array

_{T}*g*then contains the numbers of the vertices spanning the tetrahedron

^{T}*T*. During the setup step, also the vertex-to-element connectivity is required, which is simply given by

*O*(

*N*+

_{V}*N*).

_{T}*O*(

*N*) memory.

_{T}### 3.3. Sparse matrix format

### 3.4. Assembly of the system matrices

*κ*,

*μ*,

*ρ*and

*c*, the contribution of an element

*T*to the excitation and emission system is given as Therefore, we split the assembly into two steps: First, the pre-computed element matrices

*K*,

^{T}*M*and

^{T}*R*are scaled with the element’s optical parameters and the output is stored in a vector

^{T}*H*with length (4 · 4 ·

*N*). In the second step, we sum for every position in the CRS matrix the corresponding elements from

_{T}*H*. In principle, this summation can be done by multiplying a suitable sparse matrix with

*H*. However, as all non-zero elements in this matrix are 1, we implemented a specialized GPU kernel for this task, which simply sums the elements whose indices are stored in a look-up table. By using such a two-step scheme, we avoid non-local write-access and need non-local read-access only for the second summation step, which is implemented using the GPU’s texture cache. Figure 2 illustrates this procedure. As the construction of the local element matrices does only require local memory access, it can be done fully in parallel which provides almost optimal performance on both, CPU and GPU processors. The total complexity of this step is

*O*(

*N*).

_{T}### 3.5. Solution of the linear systems

*A*=

_{x}V_{x}*Q*and

*A*=

_{m}V_{m}*G*(

*c*)

*V*for the forward problem, respectively,

_{x}*A*=

_{m}W_{m}*D*and

*A*=

_{x}W_{x}*G*(

*c*)

*W*for the adjoint problem, we employ a preconditioned blocked conjugate gradient method (PBCG): instead of solving for one right hand side (i.e. a column in

_{m}*Q*) after another, we solve simultaneously for all right hand sides (blocked). The algorithm uses the PCG framework provided by the

*AGILE*library. Thus, only a new kernel had to be coded for the blocked matrix-vector-product.

*V*-cycle of a geometric multigrid preconditioner [25

25. W. L. Briggs, V. E. Henson, and S. F. McCormick, *A multigrid tutorial* (SIAM, 2000), 2nd ed. [CrossRef]

*V*-cycle has linear complexity. If the right-hand side has

*ns*columns, one application of the multigrid preconditioner requires

*O*(

*ns*·

*N*) algebraic operations.

_{V}*O*(

*ns*·

*N*) and

_{V}*O*(

*nd*·

*N*), respectively.

_{V}### 3.6. Assembly of the sensitivity

*h*is a piecewise constant function, the mass-

*M*(

*h*) and stiffness-matrices

*K*(

*h*) are just the element matrices

*M*and

^{T}*K*for the

^{T}*T*-th finite element scaled by the derivatives of the optical parameters. Since the computations for the individual elements are independent from each other, this loop can be executed in parallel.

**x**,

**y**,

**z**) = (4 × 4 × 32). The size of (4 × 4) is due to the size of the element-wise mass- and stiffness-matrices on a tetrahedral mesh. The last dimension is the number of elements which are processed in parallel and has been set to 32 to maximize the number of threads per block.

**z**instead of writing them as an array. Thus,

*n*

**is actually an array with 32 entries,**

^{z}*d*

**[**

^{z}**x**] is a 2D array with size (32 × 4) and

*k*

**(**

^{z}**x,y**) is a 3D array with size (32 × 4 × 4), for example.

*g*denotes the element-to-vertex connectivity of the

^{T}*T*-th finite element, which is cached into the vector

*d*. The algorithm caches the element’s mass- and stiffness-matrices (lines 4–6) and the forward and ajoint solutions (lines 9–10 and 12–13). Line 15 is an abbreviation for the computation of the sensitivity contribution according to Eq. (9). These contributions are summed in parallel (line 16) and the result is stored in

*S*.

*O*(

*ns*·

*nd*·

*N*) operations and the same amount of memory. Hence, this is the most time consuming and memory intensive step in the current algorithm. We would also like to note that the assembly of the sensitivity matrices

_{T}*S*can be avoided completely if it is applied on-the-fly when computing the matrix-vector products; cf [19

_{k}19. H. Egger and M. Schlottbom, “Efficient reliable image reconstruction schemes for diffuse optical tomography,” Inv. Probl. Sci. Eng. **19**, 155–180 (2011). [CrossRef]

### 3.7. Solution of the Gauß-Newton system

*N*×

_{T}*N*matrix, which would be expensive to compute and almost impossible to store. Since by construction

_{T}*S*and

_{k}*O*(

*ns*·

*nd*·

*N*) algebraic operations, whereas a multiplication with the pre-multiplied matrix

_{T}*O*(

*N*

_{T2}).

*LAPACK-BLAS*[26

26. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, *LAPACK Users’ Guide* (SIAM, 1999), 3rd ed. [CrossRef]

*AGILE*for the GPU version.

## 4. Results

### 4.1. Model problem

27. B. Dogdas, D. Stout, A. F. Chatziiouannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. **52**, 577–587 (2007). [CrossRef] [PubMed]

*ns*= 24 sources and

*nd*= 24 detectors, which were arranged on three rings around the homogeneous torso. Two spherical fluorophore inclusions with a diameter of 5 mm were placed at the height of the middle optode ring inside the body, and the aim of our test was to reconstruct these inclusions. The setup is depicted in Fig. 3.

28. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. **50**, 4225–4241 (2005). [CrossRef] [PubMed]

31. M. L. Landsman, G. Kwant, G. A. Mook, and W. G. Zijlstra, “Light-absorbing properties, stability, and spectral stabilization of indocyanine green,” J. Appl. Physiol. **40**, 575–583 (1976). [PubMed]

32. J. Schöberl, “NETGEN - an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Vis. Sci. **1**, 41–52 (1997). [CrossRef]

### 4.2. Choice of parameters

*α*and the stopping index maxit in the Gauß-Newton algorithm were chosen as follows: We conducted a series of numerical tests for ten similar phantoms, and determined a set of parameters

_{k}*α*that worked reasonably well for all samples. These parameters were then used for the actual reconstruction. For the results presented here, a geometric sequence

_{k}*α*=

_{k}*qα*

_{0}with

*α*

_{0}= 1 and

*q*= 0.2 was used. The Gauß-Newton iteration was stopped when the regularization parameter reached the pre-defined threshold

*α*= 1 × 10

_{min}^{−5}, which corresponds to maxit=8 Newton iterations.

### 4.3. Numerical results

### 4.4. Performance analysis

### 4.5. Accuracy

*L*

^{2}norm

*c*,

_{CS}*c*and

_{CD}*c*, these relative errors are The resultant images differ by less than 1% which we consider to be sufficient for practical applications.

_{GS}### 4.6. Memory requirements

*O*(

*ns*·

*nd*·

*N*)) can be reduced, i.e. it has been demonstrated in [19

_{T}19. H. Egger and M. Schlottbom, “Efficient reliable image reconstruction schemes for diffuse optical tomography,” Inv. Probl. Sci. Eng. **19**, 155–180 (2011). [CrossRef]

*O*(

*ns*·

*N*) operations. Such an approach, however, requires the solution of the adjoint problems for every multiplication with the sensitivity matrix. Comparing with the computation times in Table 2, this would drastically slow down the computations in practice, at least for the values of

_{V}*ns*,

*nd*, and

*N*,

_{T}*N*used in our test study. If the number of sources or detectors becomes very large (e.g. if a camera is used as detector), such an alternative might become favorable.

_{V}## 5. Discussion

35. D. Göddeke, H. Wobker, R. Strzodka, J. Mohd-Yusof, P. S. McCormick, and S. Turek, “Co-processor acceleration of an unmodified parallel solid mechanics code with FEASTGPU,” Int. J. Comput. Sci. Eng. **4**, 254–269 (2009). [CrossRef]

36. M. Köster, D. Göddeke, H. Wobker, and S. Turek, “How to gain speedups of 1000 on single processors with fast FEM solvers – benchmarking numerical and computational efficiency,” Tech. rep., Fakultät für Mathematik, TU Dortmund (2008). Ergebnisberichte des Instituts für Angewandte Mathematik, Nummer 382.

4. B. Zhang, X. Yang, F. Yang, X. Yang, C. Qin, D. Han, X. Ma, K. Liu, and J. Tian, “The CUBLAS and CULA based GPU acceleration of adaptive finite element framework for bioluminescence tomography,” Opt. Express **18**, 20201–20214 (2010). [CrossRef] [PubMed]

5. J. Prakash, V. Chandrasekharan, V. Upendra, and P. K. Yalavarthy, “Accelerating frequency-domain diffuse optical tomographic image reconstruction using graphics processing units,” J. Biomed. Opt. **15**, 066009 (2010). [CrossRef]

*coalesced*memory access) [33]. Therefore, the photon field matrices

*V*,

_{x}*V*and

_{m}*W*,

_{x}*W*, but also the associated source and detector matrices

_{m}*Q*and

*D*are stored in row-major order. As one column of these matrices belongs to one photon field, this means that the elements of a given photon field are interleaved. In combination with the blocked conjugate gradient method, internal performance measurements have shown that this memory layout provides speed-ups by a factor of roughly 2–3 on CPUs and up to 10 on GPUs.

*α*crucially determines the performance of the iterative reconstruction algorithm. Following the analysis of the iteratively regularized Gauß-Newton method [12, 15

_{k}15. B. Kaltenbacher, “Some Newton-type methods for the regularization of nonlinear ill-posed problems,” Inv. Probl. **13**, 729–753 (1997). [CrossRef]

*α*should decay exponentially, e.g.

_{k}*α*=

_{k}*α*

_{0}

*q*for some 0 <

^{k}*q*< 1. The choice of

*α*

_{0}and

*q*, however, depends on the problem. While in principle, appropriate values can be found by a careful analysis of the problem under investigation, a more practical way to tune these parameters is to perform a sequence of reconstructions for known solutions and select regularization parameter that work well for all cases in this training set. The same set of parameters can then be used for similar problems.

*α*

_{0}and

*q*have been fixed, and

*α*=

_{k}*α*

_{0}

*q*, then this amounts to choosing a minimal regularization parameter

^{k}*α*:=

_{min}*α*

_{maxit}. A rigorous way to choose

*α*(respectively maxit) is provided by Morozov’s discrepancy principle [37, 38], i.e. the reconstruction algorithm is stopped, as soon as ||ℳ(

_{min}*c*) – ℳ

_{k}*|| ∼*

^{δ}*δ*, where

*δ*is a measure for the measurement and model errors. Another rigorous way would be to choose

*α*∼

*δ*. Both choices however require the knowledge of the level

*δ*of the perturbations, which is usually not known in practice. A heuristic way to choose

*α*without knowledge of

_{min}*δ*is provided by the L-curve method [38,39

39. P. C. Hansen, *Rank-Defficient and Discrete Ill-Posed Problems* (SIAM, 1998). [CrossRef]

*α*) can again be obtained by “training” the algorithm on a small data set with known solutions. This strategy is probably the most practical one and also the method we utilize in our experiments.

_{min}## Acknowledgments

## References and links

1. | S. Mordon, V. Maunoury, J. M. Devoisselle, Y. Abbas, and D. Coustaud, “Characterization of tumorous and normal tissue using a pH-sensitive fluorescence indicator (5,6-carboxyfluorescein) in vivo,” J. Photochem. Photobiol. B |

2. | I. Gannot, I. Ron, F. Hekmat, V. Chernomordik, and A. Gandjbakhche, “Functional optical detection based on pH dependent fluorescence lifetime,” Lasers Surg. Med. |

3. | E. Shives, Y. Xu, and H. Jiang, “Fluorescence lifetime tomography of turbid media based on an oxygen-sensitive dye,” Opt. Express |

4. | B. Zhang, X. Yang, F. Yang, X. Yang, C. Qin, D. Han, X. Ma, K. Liu, and J. Tian, “The CUBLAS and CULA based GPU acceleration of adaptive finite element framework for bioluminescence tomography,” Opt. Express |

5. | J. Prakash, V. Chandrasekharan, V. Upendra, and P. K. Yalavarthy, “Accelerating frequency-domain diffuse optical tomographic image reconstruction using graphics processing units,” J. Biomed. Opt. |

6. | F. Knoll, M. Freiberger, K. Bredies, and R. Stollberger, “AGILE: An open source library for image reconstruction using graphics card hardware acceleration,” in “Proceedings of the 19th Scientific Meeting and Exhibition of ISMRM, Montreal, CA,” (2011), p. 2554. |

7. | S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl. |

8. | A. Joshi, W. Bangerth, and W. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express |

9. | S. R. Arridge and M. Schweiger, “The use of multiple data types in time-resolved optical absorption and scattering tomography (TOAST),” in “Mathematical Methods in Medical Imaging II,”, D. C. Wilson and J. N. Wilson, eds. (1993), |

10. | E. M. Sevick and B. Chance, “Photon migration in a model of the head measured using time and frequency domain techniques: potentials of spectroscopy and imaging,” in “Time-Resolved Spectroscopy and Imaging of Tissues,”, B. Chance and A. Katzir, eds. (1991), Proc. SPIE |

11. | S. R. Arridge, “Photon-measurement density functions. part I: Analytical forms,” Appl. Opt. |

12. | A. Bakushinsky, “The problem of the convergence of the iteratively regularized Gauss-Newton method,” Comput. Math. Math. Phys. |

13. | M. Hanke, “Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems,” Numer. Func. Anal. Optim. |

14. | D. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. |

15. | B. Kaltenbacher, “Some Newton-type methods for the regularization of nonlinear ill-posed problems,” Inv. Probl. |

16. | H. Jiang, “Frequency-domain fluorescent diffusion tomography: A finite-element-based algorithm and simulations,” Appl. Opt. |

17. | R. Roy and E. M. Sevick-Muraca, “Truncated Newtons optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation,” Opt. Express |

18. | M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss–Newton method of image reconstruction in diffuse optical tomography,” Phys. Med. Biol. |

19. | H. Egger and M. Schlottbom, “Efficient reliable image reconstruction schemes for diffuse optical tomography,” Inv. Probl. Sci. Eng. |

20. | A. Godavarty, E. M. Sevick-Muraca, and M. J. Eppstein, “Three-dimensional fluorescence lifetime tomography,” Med. Phys. |

21. | M. Freiberger, H. Egger, and H. Scharfetter, “Nonlinear inversion schemes for fluorescence optical tomography,” IEEE Trans. Biomed. Eng. |

22. | D. Göddeke, C. Becker, and S. Turek, “Integrating GPUs as fast co–processors into the parallel FE package FEAST,” in “19th Symposium Simulations Technique,”, M. Becker and H. Szczerbicka, eds., |

23. | M. M. Baskaran and R. Bordawekar, “Optimizing sparse matrix-vector multiplication on GPUs,” IBM Technical Report RC24704, IBM Ltd. (2009). |

24. | M. Liebmann, “Efficient PDE solvers on modern hardware with applications in medical and technical sciences,” Ph.D. thesis (University of Graz, 2009). |

25. | W. L. Briggs, V. E. Henson, and S. F. McCormick, |

26. | E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, |

27. | B. Dogdas, D. Stout, A. F. Chatziiouannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. |

28. | G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. |

29. | M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt. |

30. | A. Joshi, “Adaptive finite element methods for fluorescence enhanced optical tomography,” Ph.D. thesis (Texas A&M University, 2005). |

31. | M. L. Landsman, G. Kwant, G. A. Mook, and W. G. Zijlstra, “Light-absorbing properties, stability, and spectral stabilization of indocyanine green,” J. Appl. Physiol. |

32. | J. Schöberl, “NETGEN - an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Vis. Sci. |

33. | NVIDIA, |

34. | N. Bell and M. Garland, “Implementing sparse matrix-vector multiplication on throughput-oriented processors,” in “SC ’09: Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis,” (ACM, 2009), pp. 1–11. |

35. | D. Göddeke, H. Wobker, R. Strzodka, J. Mohd-Yusof, P. S. McCormick, and S. Turek, “Co-processor acceleration of an unmodified parallel solid mechanics code with FEASTGPU,” Int. J. Comput. Sci. Eng. |

36. | M. Köster, D. Göddeke, H. Wobker, and S. Turek, “How to gain speedups of 1000 on single processors with fast FEM solvers – benchmarking numerical and computational efficiency,” Tech. rep., Fakultät für Mathematik, TU Dortmund (2008). Ergebnisberichte des Instituts für Angewandte Mathematik, Nummer 382. |

37. | V. A. Morozov, “On the solution of functional equations by the method of regularization,” Soviet Math. Dokl. |

38. | H. W. Engl, M. Hanke, and A. Neubauer, |

39. | P. C. Hansen, |

40. | G. Wahba, |

**OCIS Codes**

(100.3010) Image processing : Image reconstruction techniques

(100.3190) Image processing : Inverse problems

(170.6960) Medical optics and biotechnology : Tomography

(170.7050) Medical optics and biotechnology : Turbid media

(300.6280) Spectroscopy : Spectroscopy, fluorescence and luminescence

**ToC Category:**

Image Reconstruction and Inverse Problems

**History**

Original Manuscript: June 16, 2011

Revised Manuscript: August 4, 2011

Manuscript Accepted: August 16, 2011

Published: October 28, 2011

**Citation**

Manuel Freiberger, Herbert Egger, Manfred Liebmann, and Hermann Scharfetter, "High-performance image reconstruction in fluorescence tomography on desktop computers and graphics hardware," Biomed. Opt. Express **2**, 3207-3222 (2011)

http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-2-11-3207

Sort: Year | Journal | Reset

### References

- S. Mordon, V. Maunoury, J. M. Devoisselle, Y. Abbas, and D. Coustaud, “Characterization of tumorous and normal tissue using a pH-sensitive fluorescence indicator (5,6-carboxyfluorescein) in vivo,” J. Photochem. Photobiol. B13, 307–314 (1992). [CrossRef] [PubMed]
- I. Gannot, I. Ron, F. Hekmat, V. Chernomordik, and A. Gandjbakhche, “Functional optical detection based on pH dependent fluorescence lifetime,” Lasers Surg. Med.35, 342–348 (2004). [CrossRef] [PubMed]
- E. Shives, Y. Xu, and H. Jiang, “Fluorescence lifetime tomography of turbid media based on an oxygen-sensitive dye,” Opt. Express10, 1557–1562 (2002). [PubMed]
- B. Zhang, X. Yang, F. Yang, X. Yang, C. Qin, D. Han, X. Ma, K. Liu, and J. Tian, “The CUBLAS and CULA based GPU acceleration of adaptive finite element framework for bioluminescence tomography,” Opt. Express18, 20201–20214 (2010). [CrossRef] [PubMed]
- J. Prakash, V. Chandrasekharan, V. Upendra, and P. K. Yalavarthy, “Accelerating frequency-domain diffuse optical tomographic image reconstruction using graphics processing units,” J. Biomed. Opt.15, 066009 (2010). [CrossRef]
- F. Knoll, M. Freiberger, K. Bredies, and R. Stollberger, “AGILE: An open source library for image reconstruction using graphics card hardware acceleration,” in “Proceedings of the 19th Scientific Meeting and Exhibition of ISMRM, Montreal, CA,” (2011), p. 2554.
- S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl.15, R41–R93 (1999). [CrossRef]
- A. Joshi, W. Bangerth, and W. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express12, 5402–5417 (2004). [CrossRef] [PubMed]
- S. R. Arridge and M. Schweiger, “The use of multiple data types in time-resolved optical absorption and scattering tomography (TOAST),” in “Mathematical Methods in Medical Imaging II,”, D. C. Wilson and J. N. Wilson, eds. (1993), Proc. SPIE2035, pp. 218–229. [CrossRef]
- E. M. Sevick and B. Chance, “Photon migration in a model of the head measured using time and frequency domain techniques: potentials of spectroscopy and imaging,” in “Time-Resolved Spectroscopy and Imaging of Tissues,”, B. Chance and A. Katzir, eds. (1991), Proc. SPIE1431, pp. 84–96. [CrossRef]
- S. R. Arridge, “Photon-measurement density functions. part I: Analytical forms,” Appl. Opt.34, 7395–7409 (1995). [CrossRef] [PubMed]
- A. Bakushinsky, “The problem of the convergence of the iteratively regularized Gauss-Newton method,” Comput. Math. Math. Phys.32, 1353–1359 (1992).
- M. Hanke, “Regularizing properties of a truncated Newton-CG algorithm for nonlinear inverse problems,” Numer. Func. Anal. Optim.18, 971–993 (1997). [CrossRef]
- D. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math.11, 431–441 (1963). [CrossRef]
- B. Kaltenbacher, “Some Newton-type methods for the regularization of nonlinear ill-posed problems,” Inv. Probl.13, 729–753 (1997). [CrossRef]
- H. Jiang, “Frequency-domain fluorescent diffusion tomography: A finite-element-based algorithm and simulations,” Appl. Opt.37, 5337–5343 (1998). [CrossRef]
- R. Roy and E. M. Sevick-Muraca, “Truncated Newtons optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation,” Opt. Express4, 353–371 (1999). [CrossRef] [PubMed]
- M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss–Newton method of image reconstruction in diffuse optical tomography,” Phys. Med. Biol.50, 2365–2386 (2005). [CrossRef] [PubMed]
- H. Egger and M. Schlottbom, “Efficient reliable image reconstruction schemes for diffuse optical tomography,” Inv. Probl. Sci. Eng.19, 155–180 (2011). [CrossRef]
- A. Godavarty, E. M. Sevick-Muraca, and M. J. Eppstein, “Three-dimensional fluorescence lifetime tomography,” Med. Phys.32, 992–1000 (2005). [CrossRef] [PubMed]
- M. Freiberger, H. Egger, and H. Scharfetter, “Nonlinear inversion schemes for fluorescence optical tomography,” IEEE Trans. Biomed. Eng.57, 2723–2729 (2010). [CrossRef]
- D. Göddeke, C. Becker, and S. Turek, “Integrating GPUs as fast co–processors into the parallel FE package FEAST,” in “19th Symposium Simulations Technique,”, M. Becker and H. Szczerbicka, eds., Frontiers in Simulation (SCS Publishing House, 2006), pp. 277–282.
- M. M. Baskaran and R. Bordawekar, “Optimizing sparse matrix-vector multiplication on GPUs,” IBM Technical Report RC24704, IBM Ltd. (2009).
- M. Liebmann, “Efficient PDE solvers on modern hardware with applications in medical and technical sciences,” Ph.D. thesis (University of Graz, 2009).
- W. L. Briggs, V. E. Henson, and S. F. McCormick, A multigrid tutorial (SIAM, 2000), 2nd ed. [CrossRef]
- E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users’ Guide (SIAM, 1999), 3rd ed. [CrossRef]
- B. Dogdas, D. Stout, A. F. Chatziiouannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol.52, 577–587 (2007). [CrossRef] [PubMed]
- G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol.50, 4225–4241 (2005). [CrossRef] [PubMed]
- M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt.27, 1820–1824 (1988). [CrossRef] [PubMed]
- A. Joshi, “Adaptive finite element methods for fluorescence enhanced optical tomography,” Ph.D. thesis (Texas A&M University, 2005).
- M. L. Landsman, G. Kwant, G. A. Mook, and W. G. Zijlstra, “Light-absorbing properties, stability, and spectral stabilization of indocyanine green,” J. Appl. Physiol.40, 575–583 (1976). [PubMed]
- J. Schöberl, “NETGEN - an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Vis. Sci.1, 41–52 (1997). [CrossRef]
- NVIDIA, NVIDIA CUDA Programming Guide 2.0 (NVIDIA Cooperation, 2008).
- N. Bell and M. Garland, “Implementing sparse matrix-vector multiplication on throughput-oriented processors,” in “SC ’09: Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis,” (ACM, 2009), pp. 1–11.
- D. Göddeke, H. Wobker, R. Strzodka, J. Mohd-Yusof, P. S. McCormick, and S. Turek, “Co-processor acceleration of an unmodified parallel solid mechanics code with FEASTGPU,” Int. J. Comput. Sci. Eng.4, 254–269 (2009). [CrossRef]
- M. Köster, D. Göddeke, H. Wobker, and S. Turek, “How to gain speedups of 1000 on single processors with fast FEM solvers – benchmarking numerical and computational efficiency,” Tech. rep., Fakultät für Mathematik, TU Dortmund (2008). Ergebnisberichte des Instituts für Angewandte Mathematik, Nummer 382.
- V. A. Morozov, “On the solution of functional equations by the method of regularization,” Soviet Math. Dokl.7, 414–417 (1966).
- H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems (Kluwer, 1996).
- P. C. Hansen, Rank-Defficient and Discrete Ill-Posed Problems (SIAM, 1998). [CrossRef]
- G. Wahba, Spline Models for Observational Data (SIAM, 1990).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.