OSA's Digital Library

Biomedical Optics Express

Biomedical Optics Express

  • Editor: Joseph A. Izatt
  • Vol. 4, Iss. 1 — Jan. 1, 2013
  • pp: 1–14
« Show journal navigation

Accelerated image reconstruction in fluorescence molecular tomography using dimension reduction

Xu Cao, Xin Wang, Bin Zhang, Fei Liu, Jianwen Luo, and Jing Bai  »View Author Affiliations


Biomedical Optics Express, Vol. 4, Issue 1, pp. 1-14 (2013)
http://dx.doi.org/10.1364/BOE.4.000001


View Full Text Article

Acrobat PDF (2666 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

With the development of charge-coupled device (CCD) camera based non-contact fluorescence molecular tomography (FMT) imaging systems, multi projections and densely sampled fluorescent measurements used in subsequent image reconstruction can be easily obtained. However, challenges still remain in fast image reconstruction because of the large computational burden and memory requirement in the inverse problem. In this work, an accelerated image reconstruction method in FMT using principal components analysis (PCA) is presented to reduce the dimension of the inverse problem. Phantom experiments are performed to verify the feasibility of the proposed method. The results demonstrate that the proposed method can accelerate image reconstruction in FMT almost without quality degradation.

© 2012 OSA

1. Introduction

Over the past decade, fluorescence molecular tomography (FMT) has been developed into a powerful in vivo imaging method with the advantage of non invasive and non ionizing imaging [1

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

]. The emitting fluorescent photons can be captured over the surface of small animals by a non-contact full angle FMT imaging system. Using the fluorescent measurements and appropriate model of photon propagation in biological tissues, FMT can localize and quantify fluorescent markers to resolve biological processes at molecular and cellular levels. It has been successfully applied to the studies of cancer diagnosis [2

2. E. E. Graves, R. Weissleder, and V. Ntziachristos, “Fluorescence molecular imaging of small animal tumour models,” Curr. Mol. Med. 4(4), 419–430 (2004). [CrossRef] [PubMed]

, 3

3. N. C. Deliolanis, J. Dunham, T. Wurdinger, J. L. Figueiredo, T. Bakhos, and V. Ntziachristos, “In-vivo imaging of murine tumors using complete-angle projection fluorescence molecular tomography,” J. Biomed. Opt. 14(3), 030509 (2009). [CrossRef] [PubMed]

], drug discovery [4

4. M. Rudin and R. Weissleder, “Molecular imaging in drug discovery and development,” Nat. Rev. Drug Discov. 2(2), 123–131 (2003). [CrossRef] [PubMed]

, 5

5. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discovery 7(7), 591–607 (2008). [CrossRef]

], and gene expression visualization [6

6. T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. 17(5), 545–580 (2003). [CrossRef] [PubMed]

].

Image reconstruction in FMT is a model based tomographic algorithm to compute the unknown distribution of fluorescent markers inside small animals. An accurate forward model describing the photons propagation in biological tissues is necessary for image reconstruction in FMT. Although Monte Carlo method can accurately describe the diffuse and scattering properties of photons propagation in biological tissue [7

7. L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Programs Biomed. 47(2), 131–146 (1995). [CrossRef]

], the huge computation makes it impossible to be used in fast image reconstruction. Considering the computation efficiency, diffusion equation is a widely used mathematic model in FMT. Although numerical methods are available for accurately solving diffusion equation in arbitrary geometries such as the finite difference method (FDM) [8

8. D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. 36(10), 2260–2272 (1997). [CrossRef] [PubMed]

], the finite element method (FEM) [9

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

], these methods are computationally intensive. As a fast analytical approach to solve the diffusion equation in arbitrary geometries, Kirchhoff approximation (KA) [10

10. J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Phys. Rev. E 64(5), 0519172001.

, 11

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

] is employed to solve the forward model in this work.

The inverse problem of FMT is to solve an ill-posed linear equation of which the weight matrix maps unknown distributions of the fluorescent markers inside small animals onto the fluorescent measurements over the surface. The ill-posedness of the inverse problem makes the reconstruction of FMT to be delicate and unstable to noise. To overcome the ill-posedness and get high resolution and robust reconstructed images, a large measured fluorescent data set over the surface of the small animal from 360° projections is collected by a high performance charge-coupled device (CCD) based non-contact full angle FMT imaging system [12

12. N. Deliolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360° geometry projections,” Opt. Lett. 32(4), 382–384 (2007). [CrossRef] [PubMed]

]. As each measurement corresponds to a source-detector map which is a row vector of the weight matrix, a large measured fluorescent data set leads to a large weight matrix which results in an inherently large computational burden in the inverse problem of FMT. To handle this large data set problem, the concept of data compression has been employed to encode the CCD measurements with Fourier [13

13. J. Ripoll, “Hybrid Fourier-real space method for diffuse optical tomography,” Opt. Lett. 35(5), 688–690 (2010). [CrossRef] [PubMed]

] or wavelet [14

14. T. J. Rudge, V. Y. Soloviev, and S. R. Arridge, “Fast image reconstruction in fluoresence optical tomography using data compression,” Opt. Lett. 35(5), 763–765 (2010). [CrossRef] [PubMed]

, 15

15. N. Ducros, C. D. Andrea, G. Valentini, T. Rudge, S. Arridge, and A. Bassi, “Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination,” Opt. Lett. 35(21), 3676–3678 (2010). [CrossRef] [PubMed]

, 16

16. N. Ducros, A. Bassi, G. Valentini, M. Schweiger, S. Arridge, and C. D Andrea, “Multiple-view fluorescence optical tomography reconstruction using compression of experimental data,” Opt. Lett. 36(8), 1377–1379 (2011). [CrossRef] [PubMed]

] transforms to obtain a smaller scale inverse problem for the reconstruction of FMT with all CCD measurements. A matrix-free method [17

17. A. D. Zacharopoulos, P. Svenmarker, J. Axelsson, M. Schweiger, S. R. Arridge, and S. Andersson-Engels, “A matrix-free algorithm for multiple wavelength fluorescence tomography,” Opt. Express 17(5), 3042–3051 (2009). [CrossRef]

, 18

18. A. D. Zacharopoulos, A. Garofalakis, J. Ripoll, S. R. Arridge, and S. Andersson-Engels, “Development of in-vivo fluorescence imaging with the Matrix-Free method,” J. Phys. Conf. Ser. 255(1), 012006 (2010). [CrossRef]

] is proposed to avoid explicit calculation and storage of the large weight matrix by replacing it using matrix-vector products. It can effectively reduce the computational cost and memory requirements for the FMT reconstruction. An alternative approach is to down-sample the measurements. Singular value analysis (SVA) [19

19. T. Lasser and V. Ntziachristos, “Optimization of 360° projection fluorescence molecular tomography,” Med. Image Anal. 11(4), 389–399 (2007). [CrossRef] [PubMed]

, 20

20. D. Wang, X. Liu, and J. Bai, “Analysis of fast full angle fluorescence diffuse optical tomography with beam-forming illumination,” Opt. Exp. 17(24), 21376–21395 (2009). [CrossRef]

] has been used to determine the optimal sampling interval of the measured fluorescent data compromising between reconstruction quality and computational burden. However, a large dimension of the inverse linear problem still need to be solved for image reconstruction even if the down-sampled data set with the optimal sampling interval is used.

In this work, considering the correlations of source-detector maps from the same projection, principal components analysis (PCA) is used to reduce the dimension of the weight matrix by discarding the relevant components and retaining the larger variational components. Then image reconstruction in FMT can be accelerated by solving the smaller scale inverse problem rather than the original one. Since PCA can maximally retain the useful information of the weight matrix, the computational scale can be effectively reduced almost without degrading the qualities of the reconstructed images. To verify the feasibility of the proposed method, several phantom experiments are implemented. The results demonstrate that the proposed method can accelerate image reconstruction in FMT almost without quality degradation. Besides, a comparison is made between the proposed dimension reduction method based on PCA (PCA method) and the data compression method based on wavelet transformation (wavelet transformation method) [14

14. T. J. Rudge, V. Y. Soloviev, and S. R. Arridge, “Fast image reconstruction in fluoresence optical tomography using data compression,” Opt. Lett. 35(5), 763–765 (2010). [CrossRef] [PubMed]

].

2. Method

2.1. Forward model

Considering a continuous wave (CW) point excitation source at rs, the diffusion equation with the Robin-type boundary condition [21

21. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite elementmethod for the propagation of light in scatteringmedia: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef] [PubMed]

] to describe the light transportation field in tissue is
{DG(rs,r)μaG(rs,r)=δ(rrs),2D(b)G(rs,b)n+qG(rs,b)=0
(1)
where r is the position vector in the domain, b is the position vector on the boundary, μa is the absorption coefficient, D = 1/(3(μa + μs)) is the diffusion coefficient of the tissue with reduced scattering coefficient μs, n⃗ is the outward normal vector to the surface and q is a constant associated with the ratio of optical reflective index of the inner tissue to that outside the boundary. G(rs, r) is the Green’s function describing the propagation of light from point rs to r, where rs is the point located at one transport mean free path ltr = 1/μs into the medium from the illumination spot.

The fluorescent measurement ϕ(rs, rd) detected at location rd due to an illumination spot located at rs can be formulated using the Born approximation as follows [22

22. R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging 23(4), 492–500 (2004). [CrossRef] [PubMed]

]
ϕ(rs,rd)=ΘG(rs,r)x(r)G(rd,r)d3r,
(2)
where G(rs, r) and G(rd, r) are the Green’s functions of excitation and emission light respectively, Θ is a unit-less constant taking account of the unknown gain and attenuation factors of the system. x(r) is the unknown distribution of fluorescent targets.

After the imaged object is discretized, a KA solution of Eq.1 is substituted into Eq.2
ϕ(rs,rd)=Ws,dX,
(3)
where X is a column vector representing the concentration of fluorescent targets to be reconstructed, Ws,d is a row vector representing the source detector map with the element of
Ws,d(r)=ΘΔvG(rs,r)G(rd,r),
(4)
where Δv is the volume of the discrete voxel.

For the sth projection, the sub matrix equation is formed by assembling all the source detector maps and the fluorescent measurements of the sth projection as follows
ϕs=WsX,
(5)
where Ws = {Ws,d} is the sub weight matrix constituted by all the source-detector maps of the sth projection and ϕs = {ϕ(rs)} is a column vector constituted by fluorescent measurements acquired from the sth projection.

Finally, the following matrix equation of the inverse problem is formed by assembling the sub matrix equations of all the projections as follows
ϕ=WX,
(6)
where W = {Ws} is the weight matrix mapping unknown distributions of the fluorescent markers inside the small animal onto the measured fluorescent data over the surface. ϕ = {ϕs} is a column vector constituted by fluorescent measurements from all the projections.

2.2. Kirchhoff approximation

KA, also called the tangent-plane method, is an analytical method for the diffusion equation in arbitrary geometries, which calculates the light intensity at any point of the surface by summing the homogeneous incident and reflected intensity [10

10. J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Phys. Rev. E 64(5), 0519172001.

]. The KA solution of Eq.1 for discretized situation can be expressed as follows [11

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

]
GKA(rd)=g(κ|rsrd|)+p=1N[g(κ|rprd|)+2CndDg(κ|rsrd|)np]GKA(rp)npΔS(rp)
(7)
where GKA represents the total intensity given by the KA, rp is a point on the the surface, n⃗p is the outward normal vector at point rp, ΔS(rp) is the locally planar discrete area, Cnd is a constant that takes into account the refractive index mismatch between the diffusive and the non-diffusive media, κ=μa/D is the wave number, g is the infinite homogeneous Greens function
g(κ|rsrd|)=eκ|rsrd|4πD|rsrd|,
(8)
and the surface values ∂GKA(rp)/∂np can be obtained by:
GKA(rp)np=[g(Z)g(Z+4CndD)]2CndD,
(9)
where Z = (r⃗sr⃗p) · (−n⃗p) · n⃗p. Comparing with many numerical methods for solving the diffusion equation, KA is not only able to handle arbitrary geometries, but also has the advantage of computation efficiency, because it is a linear method that does not involve matrix inversion when solving diffusion equation in arbitrary geometries. Ref. [10

10. J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Phys. Rev. E 64(5), 0519172001.

] demonstrated that KA performs with an error less than 5% when the average radius of the geometry considered is R > 3κ and is more than two orders of magnitude faster than accurate numerical methods.

2.3. Dimension reduction by PCA

As a statistical technique to reduce the dimension of multi variables, PCA linearly transforms an original set of variables into a substantially smaller set of uncorrelated variables by an orthogonal projection [23

23. I. T. Jolliffe, Principal Component Analysis, 3rd ed. (Springer-Verlag, 2002).

]. These uncorrelated variables can represent most of the information in the original set of variables. The row vectors of the sub weight matrix Ws are approximately correlative since they are corresponding to the neighboring detectors with the same light source. By discarding the correlative components and retaining the larger variational components, PCA can effectively reduce the dimension of the weight matrix.

To perform PCA, each column vector of the sub weight matrix Ws can be treated as a multivariate. The covariance matrix C of the sub weight matrix Ws (′ means the transfer operation of a matrix) can be diagonalized as follows
C=PΛP,
(10)
where P is the matrix of eigenvectors of the covariance matrix C. Λ is a diagonal matrix with the elements of the eigenvalues λ1λ2λ3 ⋯ of the covariance matrix C. The principal components Ξs of the sub weight matrix Ws can be obtain by
Ξs=PWs.
(11)
Accordingly, the left side of Eq. 5 should be multiplied by P
Γs=Pϕs.
(12)
Then a transformed sub matrix equation for the sth projection is obtained
Γs=ΞsX.
(13)
After retaining the first t large principal components and leaving out the rest less significant principal components of Ξs and Γs, a dimension reduced sub matrix equation is given as follows
Γst=ΞstX.
(14)

In this work, the cumulative percent of variance (CPVt) is used to decide the number of retained principal components t[24

24. D. A. Jackson, “Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches,” Ecology 74(8), 2204–2214 (1993). [CrossRef]

]
CPVt=i=1tλiλi.
(15)
The number of retained principal components t is determined when the CPVt reaches to a preset threshold CPVTH.

Finally, the dimension reduced matrix equation of the inverse problem is formed by assembling the dimension reduced sub matrix equations of all the projections as follows
Γt=ΞtX,
(16)
where Γt={Γst} and Ξt={Ξst}.

2.4. Tikhonov regularization method

To overcome the ill-posedness, linear system Eq.16 is solved by iterated Tikhonov regularization method [25

25. M. Hanke and C. W. Groetsch, “Nonstationary iterated Tikhonov regularization,” J. Optim. Theor. Appl. 98(1), 37–53 (1998). [CrossRef]

, 26

26. V. Faber, A. Manteuffel, A. B. White Jr., and G. M. Wing, “Asymptotic behavior of singular values and functions of certain convolution operators,” Comput. Math. Appl. 12, 37–52 (1989).

]
Xn=Xn1(ΞtΞt+αI)1[Ξt(Ξt*Xn1Γt)],
(17)
Considering the row of the reduced matrix number is much smaller than column number, we using the right inverse formation as follow [27

27. C. R. Rao and S. K. Mitra, Generalized Inverse of Matrices and Its Applications (Wiley, New York, 1971).

]
Xn=Xn1Ξt[(ΞtΞt+αI)1(Ξt*Xn1Γt)],
(18)
where α is the regularization parameter. The initial X0 is set to X0 = 0, the number of iterations is set to 30, and the regularization parameter is set to 10−5trt Ξ′t) for the following phantom experiments. Note that (Ξt Ξ′t + αI)−1 is the most time-consuming process in Eq.18, which involves matrix multiplication and matrix inversion. The time cost in calculating (Ξt Ξ′t + αI)−1 is reduced significantly when the row number of original weight matrix W is reduced into Ξt. Fortunately, (Ξt Ξ′t + αI)−1 needs to be calculated only once, as it is independent on iteration parameter n. Except for the matrix inversion, the other processes only contain matrix-vector multiplications and cost little time.

3. Experimental setup and materials

3.1. Experimental setup

Phantom experiments were conducted based on the experimental setup shown in Fig. 1, which is detailed in [28

28. F. Liu, X. Liu, D. F. Wang, B. Zhang, and J. Bai, “Parallel Excitation Based Fluorescence Molecular Tomography System for Whole-Body Simultaneous Imaging of Small Animals,” Ann. Biomed. Eng. 38(11), 3440–3448 (2010). [CrossRef] [PubMed]

]. A 250W Halogen lamp (7ILT250, 7-star, Beijing, China), was employed as the excitation light source. The spot excitation source was provided by the excitation light traveling through a 775±23nm band-pass filter (FF01-775/46-25, Semrock, Rochester, NY, USA) coupled with an optical fiber. An electron multiplying charge-coupled device (EMCCD) camera (iXon DU-897, Andor Technologies, Belfast, Northern Ireland) coupled with a Nikkor 60 mm, f/2.8D lens (Nikon,Melville, NY, USA) was placed on the opposite side of the imaged object to collect photons propagating through it in the transillumination mode. The EMCCD was cooled to −70°C to reduce dark noise. An 840 ± 6nm band-pass filter (FF01-840/12-25, Semrock, Rochester, NY, USA) was placed in front of the camera for fluorescence detection. When collecting the white light images for recovering the geometry of the imaged object, an incandescent lamp was employed to serve as the excitation source. The exposure time for collecting fluorescent images with 2 × 2 binning was 2s. 36 fluorescent images were obtained in 10° steps. 72 white light images was collected every 5° during a 360° rotation.

Fig. 1 The schematic of the experimental setup.

3.2. Phantom studies

In the phantom studies, a transparent glass cylinder with a diameter of 2.4cm, height of 3.8cm filled with 1% intralipid (μa=0.02cm−1 and μs=10cm−1) [29

29. G. Q. Yu, T. Durduran, C. Zhou, H. W. Wang, M. E. Putt, H. M Saunders, C. M. Sehgal, E. Glatstein, A. G. Yodh, and T. M. Busch, “Noninvasive monitoring of murine tumor blood flow during and after photodynamic therapy provides early assessment of therapeutic efficacy,” Clin Cancer Res. 1(9), 3543–3552 (2005). [CrossRef]

] was employed as the phantom. For the single fluorescent target case, a small transparent glass tube (0.5cm in diameter) filled with 10μL ICG with the concentration of 1.3μM, was immersed in the phantom at (−0.1cm, 0cm, 1.5cm). For the double fluorescent targets case, two small transparent glass tubes (0.4cm in diameter) filled with 10μL ICG with the concentration of 6.5μM were immersed in the phantom with the edge to edge distance of 0.5cm (one was at (−0.1cm, 0.5cm, 1.6cm), the other at (0cm, −0.4cm, 1.6cm)).

4. Results and discussion

4.1. Phantom results

Figures 2 and 3 are the reconstructed results of the single and double fluorescent target cases, respectively. CPVTH = 90% is used to determine the number of retained large principal components of each sub weight matrix Ws. There are almost no visual differences in both three dimensional (3D) views and sections of the results reconstructed from the reduced weight matrix using PCA and the original weight matrix in the two experiments. The time cost and sizes of dimension reduced and original weight matrices for the two experiments are listed in Table 1. As shown in the table, both the methods need to solve forward problem with the same scale, so the time cost in the calculation of forward problem is the same. In the assembling process, due to the extra compression or dimension reduction operation, assembling by the reduce weight matrix using PCA costs more time than the original weight matrix. In the matrix inversion process, using the reduced weight matrix can save lots of time compared with the original weight matrix. As a result, the total reconstruction time of the proposed method is less than one minute, which is about 1/8 of the time cost by using the original weight matrix. These results demonstrate that the proposed method can effectively accelerate the reconstruction almost without degrading the qualities of the reconstructed images.

Fig. 2 Reconstruction results of the single fluorescent target case. (a) Section and (d) 3D view of true fluorescent target. (b) Section and (e) 3D view of reconstructed result using the original weight matrix. (c) Section and (f) 3D view of reconstructed result using the dimension reduced weight matrix. All the images are normalized by the maximum of the results.
Fig. 3 Reconstruction results of the double fluorescent targets case. (a) Section and (d) 3D view of true fluorescent target. (b) Section and (e) 3D view of reconstructed result using the original weight matrix. (c) Section and (f) 3D view of reconstructed result using the dimension reduced weight matrix. All the images are normalized by the maximum of the results.

Table 1. Quantification analysis of computational scales

table-icon
View This Table

4.2. Influence of the number of retained principal components

To further analyze the influence of the number of retained principal components of each sub weight matrix Ws on the qualities of reconstructed images and computation time, we investigated the reconstructed results using different CPVTH. The errors sqrt(||XreducedXoriginal||/||Xoriginal||) are calculated to quantify the differences between the results reconstructed from the dimension reduced weight matrix and the original weight matrix, where Xreduced and Xoriginal are the results reconstructed using the dimension reduced weight matrix and the original weight matrix, respectively.

The reconstructed results, corresponding to the errors and computation time using different CPVTH for the single and double fluorescent target cases are shown in Figs. 4 and 5, respectively. Figures 4 (a) and 5 (a) are the sections of reconstructed results using the original weight matrix. Figures 4 (b)–(h) and 5 (b)–(h) are the sections of reconstructed results obtained using the dimension reduction weight matrix when CPVTH are 98%, 95%, 90%, 85%, 80%, 75% and 70%, respectively. In the single fluorescent target case, it is hard to find visual differences in the results reconstructed from the original weight matrix and the dimension reduced weight matrices using different CPVTH. In the double fluorescent targets case, there appears some stains between the targets when CPVTH = 70% but no distinct differences for other larger CPVTH. That is because resolving a complex distribution of multi fluorescent targets usually needs more projected information. To quantify the above results, Figs. 4 (i) and 5 (i) show the errors and computation time, respectively. It can be seen that increasing CPVTH improves the reconstruction image quality and reduces the reconstruction errors at the cost of increasing computation time. A compromise method is to choose a appropriate CPVTH considering both image quality and computation time. The CPVTH of 90% is suggested as a compromise between reconstruction image quality and computation time.

Fig. 4 Influence of the CPVTH on the reconstructed images in the single fluorescent target case. (a) Section of the reconstructed result using the original weight matrix. (b)–(h) Sections of reconstructed results when the CPVTH are 98%, 95%, 90%, 85%, 80%, 75% and 70% respectively. (i) Errors between results using the original weight matrix and the dimension reduction matrix when different CPVTH are used to determine the number of retained principal components, and computation time as a function of different CPVTH. All the images are normalized by the maximum of the results.
Fig. 5 Influence of the CPVTH on the reconstructed images in the double fluorescent targets case. (a) Section of reconstructed result using the original weight matrix. (b)–(h) Sections of reconstructed results when the CPVTH are 98%, 95%, 90%, 85%, 80%, 75% and 70% respectively. (i) Errors between results using the original weight matrix and the dimension reduction matrix when different CPVTH are used to determine the number of retained principal components, and computation time as a function of different CPVTH. All the images are normalized by the maximum of the results.

4.3. Comparison with data compression method

In this section a comparison is made between the proposed PCA method and the wavelet transformation method. To make a reasonable comparison between these two methods, the quality of reconstructed image and the computational time are investigated when the reduced weight matrices with same scale are obtained with the PCA method and the wavelet transformation method, respectively. First we define compress rate (CR) which is the ratio of the row number of the original weight matrix against that of the reduced weight matrix. Then we can change the CPVTH of the PCA method and the retained wavelet coefficients of the wavelet transformation method to acquired a same CR of the original weight matrix.

The reconstructed results using the PCA method and the wavelet transformation method are shown in Figs. 6 and 7 for the single and double fluorescent target cases, respectively. With the increase of CR, the quality of reconstructed images obtained with the PCA method is degraded more slowly than that obtained with the wavelet transformation method. When the compress rate is large (CR=8.0, 12.9), both in the single and double fluorescent target cases, there are obvious distortions of the reconstructed fluorescent targets in the reconstructed images obtained from the wavelet transformation method. However, the reconstructed images obtained from the PCA method almost remain unchanged even at large compress rate. It demonstrate that the reduced weight matrix using PCA can retain more useful information than the wavelet transformation. Figures 6 (k) and 7 (k) show the errors and computation time for the PCA method and the wavelet transformation method, respectively. Using the same CR, the error of reconstructed result obtained with the PCA method is significantly smaller than that with the wavelet transformation method, although the computation time of the PCA method is slightly larger than that of the wavelet transformation method. Furthermore, the black dash-dot lines in Figs. 6 (k) and 7 (k) illustrate the comparisons of the computation time when the errors are approximate. In the single fluorescent target case, when the errors are 0.2357 for the PCA method (CR=8.0, Fig. 6 (d)) and 0.2346 for the wavelet transformation method (CR=3.5, Fig. 6 (g)), the computation time of the PCA method is 24.4s while that of the wavelet transformation method is 34.1s. In the double fluorescent targets case, when the errors are 0.4168 for the PCA method (CR=8.0, Fig. 7 (d)) and 0.4046 for the wavelet transformation method (CR=3.5, Fig. 7 (g)), the computation time of the PCA method is 25.9s while that of the wavelet transformation method is 35.8s. It means that the PCA method is faster than the wavelet transformation method when the qualities of reconstructed results are approximate.

Fig. 6 Comparison between the PCA method and the wavelet transformation method in the single fluorescent target case. (a) Section of true fluorescent target. (f) Section of the reconstructed result using the original weight matrix. (b)–(e) Sections of reconstructed results obtained from the PCA method when the CR of the reduced weight matrix is 3.5, 4.7, 8.0 and 12.9, respectively. (g)–(j) Sections of reconstructed results obtained from the wavelet transformation method when the CR of the reduced weight matrix is 3.5, 4.7, 8.0 and 12.9, respectively. (k) Errors and computation time as a function of different CR for the PCA method and the wavelet transformation method. All the images are normalized by the maximum of the results.
Fig. 7 Comparison between the PCA method and the wavelet transformation method in the double fluorescent targets case. (a) Section of true fluorescent target. (f) Section of the reconstructed result using the original weight matrix. (b)–(e) Sections of reconstructed results obtained from the PCA method when the CR of the reduced weight matrix is 3.5, 4.7, 8.0 and 12.9, respectively. (g)–(j) Sections of reconstructed results obtained from the wavelet transformation method when the CR of the reduced weight matrix is 3.5, 4.7, 8.0 and 12.9, respectively. (k) Errors and computation time as a function of different CR for the PCA method and the wavelet transformation method. All the images are normalized by the maximum of the results.

5. Conclusion

In this study, an accelerated image reconstruction method in FMT based on reducing the dimension of the weight matrix using PCA is presented. It can effectively reduce the computation scale and maximally retain useful information of the original weight matrix. The results of phantom experiments demonstrate that the proposed method can accelerate FMT reconstruction almost without quality degradation.

Due to the ill-posedness of the inverse problem in FMT, a large measured fluorescent data set is helpful to reconstruct a high quality image especially for high resolution situations [30

30. V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A 19(3), 558–566 (2002). [CrossRef]

]. But the large measured data set brings in a computational burden for the forward problem and inverse problem of FMT. Since these source-detector maps for one projection are generated based on the same excitation source and neighboring detectors, there are correlations among these vectors in each sub weight matrix. So PCA is a desired approach to reduce the computational scale and retain most information by reducing the dimension of each sub weight matrix. When an appropriate number of principal components of each sub weight matrix are retained, there is almost no quality degradation of images reconstructed from dimensional reduced weight matrix but can significantly save calculation time. The comparison between the proposed PCA method and the wavelet transformation method demonstrates that PCA can retain more useful information than the wavelet transformation.

It is worth noting that the reconstruction will be accelerated more effectively with the proposed method when a larger measured fluorescent data set is used in FMT reconstruction, which will emerge when measurements are acquired from more projections with a small sampling interval. FMT reconstruction with all the CCD measurements without sampling can improve the qualities of reconstructed images. Due to the correlations in the large weight matrix, the proposed method can also be extended to solve the reconstruction when all the CCD measurements are used.

Acknowledgments

This work is supported by the National Basic Research Program of China (973) under Grant No. 2011CB707701; the National Natural Science Foundation of China under Grant No. 81071191, 81271617, 81201160, 61271131; the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114; National Science and technology support program under Grant No. 2012BAI23B00; Tsinghua University Initiative Scientific Research Program; China’s Post-doctoral Science Foundation Grant No. 2012M510463.

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,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

2.

E. E. Graves, R. Weissleder, and V. Ntziachristos, “Fluorescence molecular imaging of small animal tumour models,” Curr. Mol. Med. 4(4), 419–430 (2004). [CrossRef] [PubMed]

3.

N. C. Deliolanis, J. Dunham, T. Wurdinger, J. L. Figueiredo, T. Bakhos, and V. Ntziachristos, “In-vivo imaging of murine tumors using complete-angle projection fluorescence molecular tomography,” J. Biomed. Opt. 14(3), 030509 (2009). [CrossRef] [PubMed]

4.

M. Rudin and R. Weissleder, “Molecular imaging in drug discovery and development,” Nat. Rev. Drug Discov. 2(2), 123–131 (2003). [CrossRef] [PubMed]

5.

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discovery 7(7), 591–607 (2008). [CrossRef]

6.

T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. 17(5), 545–580 (2003). [CrossRef] [PubMed]

7.

L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Programs Biomed. 47(2), 131–146 (1995). [CrossRef]

8.

D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. 36(10), 2260–2272 (1997). [CrossRef] [PubMed]

9.

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

10.

J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Phys. Rev. E 64(5), 0519172001.

11.

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

12.

N. Deliolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360° geometry projections,” Opt. Lett. 32(4), 382–384 (2007). [CrossRef] [PubMed]

13.

J. Ripoll, “Hybrid Fourier-real space method for diffuse optical tomography,” Opt. Lett. 35(5), 688–690 (2010). [CrossRef] [PubMed]

14.

T. J. Rudge, V. Y. Soloviev, and S. R. Arridge, “Fast image reconstruction in fluoresence optical tomography using data compression,” Opt. Lett. 35(5), 763–765 (2010). [CrossRef] [PubMed]

15.

N. Ducros, C. D. Andrea, G. Valentini, T. Rudge, S. Arridge, and A. Bassi, “Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination,” Opt. Lett. 35(21), 3676–3678 (2010). [CrossRef] [PubMed]

16.

N. Ducros, A. Bassi, G. Valentini, M. Schweiger, S. Arridge, and C. D Andrea, “Multiple-view fluorescence optical tomography reconstruction using compression of experimental data,” Opt. Lett. 36(8), 1377–1379 (2011). [CrossRef] [PubMed]

17.

A. D. Zacharopoulos, P. Svenmarker, J. Axelsson, M. Schweiger, S. R. Arridge, and S. Andersson-Engels, “A matrix-free algorithm for multiple wavelength fluorescence tomography,” Opt. Express 17(5), 3042–3051 (2009). [CrossRef]

18.

A. D. Zacharopoulos, A. Garofalakis, J. Ripoll, S. R. Arridge, and S. Andersson-Engels, “Development of in-vivo fluorescence imaging with the Matrix-Free method,” J. Phys. Conf. Ser. 255(1), 012006 (2010). [CrossRef]

19.

T. Lasser and V. Ntziachristos, “Optimization of 360° projection fluorescence molecular tomography,” Med. Image Anal. 11(4), 389–399 (2007). [CrossRef] [PubMed]

20.

D. Wang, X. Liu, and J. Bai, “Analysis of fast full angle fluorescence diffuse optical tomography with beam-forming illumination,” Opt. Exp. 17(24), 21376–21395 (2009). [CrossRef]

21.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite elementmethod for the propagation of light in scatteringmedia: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef] [PubMed]

22.

R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging 23(4), 492–500 (2004). [CrossRef] [PubMed]

23.

I. T. Jolliffe, Principal Component Analysis, 3rd ed. (Springer-Verlag, 2002).

24.

D. A. Jackson, “Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches,” Ecology 74(8), 2204–2214 (1993). [CrossRef]

25.

M. Hanke and C. W. Groetsch, “Nonstationary iterated Tikhonov regularization,” J. Optim. Theor. Appl. 98(1), 37–53 (1998). [CrossRef]

26.

V. Faber, A. Manteuffel, A. B. White Jr., and G. M. Wing, “Asymptotic behavior of singular values and functions of certain convolution operators,” Comput. Math. Appl. 12, 37–52 (1989).

27.

C. R. Rao and S. K. Mitra, Generalized Inverse of Matrices and Its Applications (Wiley, New York, 1971).

28.

F. Liu, X. Liu, D. F. Wang, B. Zhang, and J. Bai, “Parallel Excitation Based Fluorescence Molecular Tomography System for Whole-Body Simultaneous Imaging of Small Animals,” Ann. Biomed. Eng. 38(11), 3440–3448 (2010). [CrossRef] [PubMed]

29.

G. Q. Yu, T. Durduran, C. Zhou, H. W. Wang, M. E. Putt, H. M Saunders, C. M. Sehgal, E. Glatstein, A. G. Yodh, and T. M. Busch, “Noninvasive monitoring of murine tumor blood flow during and after photodynamic therapy provides early assessment of therapeutic efficacy,” Clin Cancer Res. 1(9), 3543–3552 (2005). [CrossRef]

30.

V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A 19(3), 558–566 (2002). [CrossRef]

OCIS Codes
(100.3190) Image processing : Inverse problems
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(170.3880) Medical optics and biotechnology : Medical and biological imaging
(170.6960) Medical optics and biotechnology : Tomography
(290.1990) Scattering : Diffusion
(290.7050) Scattering : Turbid media

ToC Category:
Image Reconstruction and Inverse Problems

History
Original Manuscript: September 12, 2012
Revised Manuscript: November 20, 2012
Manuscript Accepted: November 27, 2012
Published: December 5, 2012

Citation
Xu Cao, Xin Wang, Bin Zhang, Fei Liu, Jianwen Luo, and Jing Bai, "Accelerated image reconstruction in fluorescence molecular tomography using dimension reduction," Biomed. Opt. Express 4, 1-14 (2013)
http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-4-1-1


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol.23(3), 313–320 (2005). [CrossRef] [PubMed]
  2. E. E. Graves, R. Weissleder, and V. Ntziachristos, “Fluorescence molecular imaging of small animal tumour models,” Curr. Mol. Med.4(4), 419–430 (2004). [CrossRef] [PubMed]
  3. N. C. Deliolanis, J. Dunham, T. Wurdinger, J. L. Figueiredo, T. Bakhos, and V. Ntziachristos, “In-vivo imaging of murine tumors using complete-angle projection fluorescence molecular tomography,” J. Biomed. Opt.14(3), 030509 (2009). [CrossRef] [PubMed]
  4. M. Rudin and R. Weissleder, “Molecular imaging in drug discovery and development,” Nat. Rev. Drug Discov.2(2), 123–131 (2003). [CrossRef] [PubMed]
  5. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discovery7(7), 591–607 (2008). [CrossRef]
  6. T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev.17(5), 545–580 (2003). [CrossRef] [PubMed]
  7. L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Programs Biomed.47(2), 131–146 (1995). [CrossRef]
  8. D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt.36(10), 2260–2272 (1997). [CrossRef] [PubMed]
  9. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys.20(2), 299–309 (1993). [CrossRef] [PubMed]
  10. J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Phys. Rev. E64(5), 0519172001.
  11. J. Ripoll, M. Nieto-Vesperinas, R. Weissleder, and V. Ntziachristos, “Fast analytical approximation for arbitrary geometries in diffuse optical tomography,” Opt. Lett.27(7), 527–529 (2002). [CrossRef]
  12. N. Deliolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360° geometry projections,” Opt. Lett.32(4), 382–384 (2007). [CrossRef] [PubMed]
  13. J. Ripoll, “Hybrid Fourier-real space method for diffuse optical tomography,” Opt. Lett.35(5), 688–690 (2010). [CrossRef] [PubMed]
  14. T. J. Rudge, V. Y. Soloviev, and S. R. Arridge, “Fast image reconstruction in fluoresence optical tomography using data compression,” Opt. Lett.35(5), 763–765 (2010). [CrossRef] [PubMed]
  15. N. Ducros, C. D. Andrea, G. Valentini, T. Rudge, S. Arridge, and A. Bassi, “Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination,” Opt. Lett.35(21), 3676–3678 (2010). [CrossRef] [PubMed]
  16. N. Ducros, A. Bassi, G. Valentini, M. Schweiger, S. Arridge, and C. D Andrea, “Multiple-view fluorescence optical tomography reconstruction using compression of experimental data,” Opt. Lett.36(8), 1377–1379 (2011). [CrossRef] [PubMed]
  17. A. D. Zacharopoulos, P. Svenmarker, J. Axelsson, M. Schweiger, S. R. Arridge, and S. Andersson-Engels, “A matrix-free algorithm for multiple wavelength fluorescence tomography,” Opt. Express17(5), 3042–3051 (2009). [CrossRef]
  18. A. D. Zacharopoulos, A. Garofalakis, J. Ripoll, S. R. Arridge, and S. Andersson-Engels, “Development of in-vivo fluorescence imaging with the Matrix-Free method,” J. Phys. Conf. Ser.255(1), 012006 (2010). [CrossRef]
  19. T. Lasser and V. Ntziachristos, “Optimization of 360° projection fluorescence molecular tomography,” Med. Image Anal.11(4), 389–399 (2007). [CrossRef] [PubMed]
  20. D. Wang, X. Liu, and J. Bai, “Analysis of fast full angle fluorescence diffuse optical tomography with beam-forming illumination,” Opt. Exp.17(24), 21376–21395 (2009). [CrossRef]
  21. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite elementmethod for the propagation of light in scatteringmedia: Boundary and source conditions,” Med. Phys.22(11), 1779–1792 (1995). [CrossRef] [PubMed]
  22. R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging23(4), 492–500 (2004). [CrossRef] [PubMed]
  23. I. T. Jolliffe, Principal Component Analysis, 3rd ed. (Springer-Verlag, 2002).
  24. D. A. Jackson, “Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches,” Ecology74(8), 2204–2214 (1993). [CrossRef]
  25. M. Hanke and C. W. Groetsch, “Nonstationary iterated Tikhonov regularization,” J. Optim. Theor. Appl.98(1), 37–53 (1998). [CrossRef]
  26. V. Faber, A. Manteuffel, A. B. White, and G. M. Wing, “Asymptotic behavior of singular values and functions of certain convolution operators,” Comput. Math. Appl.12, 37–52 (1989).
  27. C. R. Rao and S. K. Mitra, Generalized Inverse of Matrices and Its Applications (Wiley, New York, 1971).
  28. F. Liu, X. Liu, D. F. Wang, B. Zhang, and J. Bai, “Parallel Excitation Based Fluorescence Molecular Tomography System for Whole-Body Simultaneous Imaging of Small Animals,” Ann. Biomed. Eng.38(11), 3440–3448 (2010). [CrossRef] [PubMed]
  29. G. Q. Yu, T. Durduran, C. Zhou, H. W. Wang, M. E. Putt, H. M Saunders, C. M. Sehgal, E. Glatstein, A. G. Yodh, and T. M. Busch, “Noninvasive monitoring of murine tumor blood flow during and after photodynamic therapy provides early assessment of therapeutic efficacy,” Clin Cancer Res.1(9), 3543–3552 (2005). [CrossRef]
  30. V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A19(3), 558–566 (2002). [CrossRef]

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.


Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited