An adaptive Tikhonov regularization method for ... - Springer Link
Med Biol Eng Comput (2013) 51:849–858 DOI 10.1007/s11517-013-1054-5
ORIGINAL ARTICLE
An adaptive Tikhonov regularization method for fluorescence molecular tomography Xu Cao • Bin Zhang • Xin Wang • Fei Liu Ke Liu • Jianwen Luo • Jing Bai
•
Received: 8 November 2012 / Accepted: 23 February 2013 / Published online: 16 March 2013 Ó International Federation for Medical and Biological Engineering 2013
Abstract The high degree of absorption and scattering of photons propagating through biological tissues makes fluorescence molecular tomography (FMT) reconstruction a severe ill-posed problem and the reconstructed result is susceptible to noise in the measurements. To obtain a reasonable solution, Tikhonov regularization (TR) is generally employed to solve the inverse problem of FMT. However, with a fixed regularization parameter, the Tikhonov solutions suffer from low resolution. In this work, an adaptive Tikhonov regularization (ATR) method is presented. Considering that large regularization parameters can smoothen the solution with low spatial resolution, while small regularization parameters can sharpen the solution with high level of noise, the ATR method adaptively updates the spatially varying regularization parameters during the iteration process and uses them to penalize the solutions. The ATR method can adequately sharpen the feasible region with fluorescent probes and smoothen the region without fluorescent probes resorting to no complementary priori information. Phantom experiments are performed to verify the feasibility of the proposed method. The results demonstrate that the proposed method can X. Cao B. Zhang X. Wang F. Liu J. Bai (&) Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China e-mail: [email protected]; [email protected] K. Liu Division of Electronics and Information Technology, National Institute of Metrology, Beijing 100013, People’s Republic of China J. Luo Center for Biomedical Imaging Research and Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China
improve the spatial resolution and reduce the noise of FMT reconstruction at the same time. Keywords Adaptive regularization Fluorescence molecular tomography regularization Reconstruction Tikhonov
1 Introduction Fluorescence molecular tomography (FMT) can reveal biological processes in vivo by three-dimensional (3D) visualization of specific fluorescent probes [16]. With the advantages of low cost, non-invasive and non-ionizing, FMT has been developed into a powerful in vivo preclinical imaging method for studies of cancer diagnosis [6, 14], drug discovery [21], and gene expression visualization [13]. As a typical inverse problem, the high degree of absorption and scattering of photons propagating through biological tissues makes the FMT reconstruction a severe ill-posed problem. It means that the solution is susceptible to noise in the measurements. Tikhonov regularization (TR) is a widespread method to deal with the ill-posed inverse problems, in which an L2-norm constraint term is added into the data-fitting term to improve the stability [26]. Regularization parameter is important for regularization method in controlling the smooth degree of the reconstructed image. Generally, a large regularization parameter results in a smooth image with low spatial resolution, while a small regularization parameter results in high spatial resolution image with lots of high-frequency noise. To make a compromise between the spatial resolution and signal-to-noise ratio of the reconstructed image, several methods are developed to choose an optimal
123
850
regularization parameter, such as the L-curve method [8, 7], the cross-validation method [5], and the discrepancy principle [22]. However, in practice, the selected regularization parameter usually leads to a poor spatial resolution in the reconstructed FMT images, and to make matters worse, the selected regularization parameter is not the optimal one due to the serious ill-posed nature and the model mismatch of FMT. So the optimal regularization parameter is often empirically chosen in practical applications [4, 3, 20]. To improve the quality of the reconstructed image using single constant regularization parameter, several adaptive regularization methods are proposed to improve the resolution and contrast of the reconstructed image in diffuse optical tomography (DOT) [9, 11, 15, 17, 25]. Similarly, the use of spatially varying regularization parameters according to the additional information obtained from multi-spectral fluorescent emission recordings [1] or the prior anatomical information obtained from CT images [10] in the reconstruction has been reported to improve image quality in FMT. But these two methods need complementary information which increases the complexity of imaging system. Due to the differences between FMT and DOT, the adaptive regularization methods proposed for DOT reconstruction cannot be applied in FMT reconstruction, and to the best of our knowledge, there is not an appropriate adaptive regularization method independent of complementary information for FMT reconstruction. In this study, we modify the conventional TR and propose an adaptive Tikhonov regularization (ATR) without using the complementary information, which dynamically updates the spatially varying regularization parameters based on results of previous iterations. The principle of this method is using different regularization parameters to spatially penalize different voxels of the reconstructed image. Large regularization parameters are used to spatially penalize the voxels which have no fluorescent probes and small regularization parameters are used to penalize the voxels which have fluorescent probes. As a result, the region with fluorescent probes can be adequately sharpened and the region without fluorescent probes can be adequately smoothened. To adaptively update the regularization parameters, the iteration starts with a large regularization parameter which is much larger than the empirically chosen one. During iteration processes, result of each iteration is used as the prior information to decide the changes of the regularization parameters in the next iteration. The regularization parameters corresponding to the zero value voxels stay the same and those corresponding to the non-zero value voxels are reduced. So the region with fluorescent probes is sharpened by smaller and smaller regularization parameters and the region without fluorescent target distribution is smoothened by large regularization parameters. Finally a
123
Med Biol Eng Comput (2013) 51:849–858
robust FMT reconstruction with high spatial resolution can be acquired. To examine the improvements of spatial resolution and reduction of the noise in FMT reconstruction, phantom experiments are implemented, respectively. Furthermore, the influences of important parameters of the proposed method on the reconstructed image are discussed.
2 Methods 2.1 Forward model Considering a continuous wave (CW) point excitation source at rs, the diffusion equation with the Robin-type boundary condition [24] to describe the light transportation field in tissue is r DrGðrs ; rÞ la Gðrs ; rÞ ¼ dðr rs Þ;
ð1Þ
where r is the position vector in the domain, b is the position vector on the boundary, la is the absorption coefficient, D = 1/(3(l a ? l0 s)) is the diffusion coefficient of the tissue with the reduced scattering coefficient l0 s. 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/l0 s into the medium from the illumination spot. The Robin boundary is given by 2DðbÞ
oGðrs ; bÞ þ qGðrs ; bÞ ¼ 0; on
ð2Þ
where 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. The fluorescent measurements /(rs, rd) detected at location rd due to an illumination spot located at rs can be formulated as follows, using the first Rytov approximation [23] Z /ðrs ; rd Þ ¼ H Gðrs ; rÞxðrÞGðrd ; rÞd3 r; ð3Þ where G(rs, r) and G(rd, r) are the Green’s functions of excitation and emission light respectively, H is a unit-less constant taking account of the unknown gain and attenuation factors of the system. x(r) is the distribution of fluorescent probe representing the fluorochrome concentration at location r multiplied by the fluorescent yield. After the imaged object is discretized, a Kirchhoff Approximation (KA) [19, 18] solution of Eq. (1) is substituted into Eq. (3) /ðrs ; rd Þ ¼ ws;d X;
ð4Þ
where the column vector X is the discrete formation of x(r) representing the distribution of fluorescent probe to be
Med Biol Eng Comput (2013) 51:849–858
851
reconstructed, ws,d is the source detector map with the element of ws;d ðrÞ ¼ HDvGðrs ; rÞGðrd ; rÞ;
pffi computed by c ¼ ½Nastop and N is the number of iterations. The reconstruction strategies can be summarized as follows
ð5Þ
where Dv is the volume of the discrete voxel. Finally, the following matrix equation of the inverse problem is formed by assembling all the source detector maps as follows U ¼ WX;
ð6Þ
where W = {ws,d} is the weight matrix mapping unknown distributions of the fluorescent probes inside the small animal onto the measured fluorescent data over the surface. U ¼ f/ðrs ; rd Þg is a column vector constituted by fluorescent measurements.
2.3 Experimental setup
where K is a diagonal matrix with diagonal elements of an(r) which is updated according to the previous iteration result Xn-1. The negative elements of Xn-1 are physically impossible, so those are reset to zero. an1 ðrÞ Xn1 ðrÞ ¼ 0 an ðrÞ ¼ ð9Þ can1 ðrÞ Xn1 ðrÞ [ 0
Phantom experiments were conducted based on the FMT system previously developed in our laboratory [12]. The schematic diagram of the system is shown in Fig. 1. The imaged object was fixed onto the rotation stage, which can rotate around its longitudinal axis for full-angle projection acquisitions. An appropriate rotation speed (6°/s) was selected to avoid skew and internal organs movement during the in vivo imaging process. The point excitation source was provided by a 250 W Halogen lamp (7ILT250, 7-star, Beijing, China). The excitation light traveled through a 775 ± 23 nm 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 ± 6 nm 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
where 0 \ c \ 1 is a shrinkage factor of the regularization parameter. In this equation, an(r) is spatially varying and adaptively updated based on the previous iteration result Xn-1. an(r) corresponding to zero value elements (i.e., voxels have no fluorescent probes) of Xn-1 stay the same, while an(r) corresponding to the non-zero value ones (i.e., voxels have fluorescent probes) are reduced. In the next iteration, a larger regularization parameter is used to spatially penalize the voxels which have no probes and a small one is used to spatially penalize the voxels which have probes. The initial regularization parameter is set as a0 = 1. So the first solution is equivalent to Eq. (7). During the iterations, the regularization parameters corresponding to the voxels with fluorescent probes get smaller and smaller. To avoid unstable solutions caused by too small regularization parameters, stop regularization parameter is set as astop. So the shrinkage factor can be
Fig. 1 The schematic of the experimental setup
2.2 Adaptive Tikhonov regularization method The matrix W on the right of Eq. (6) is generally ill-posed and TR is used to obtain a stable solution X ¼ ½W T W þ trðW T WÞaI1 W T U;
ð7Þ
where I is the identity matrix, tr(*) is trace of a matrix and a is the regularization parameter to stabilize the solution. In this study, Eq. (7) is modified into an iteration process using an adaptive spatially varying regularization parameter an(r) rather than a single scalar value a Xn ¼ ½W T W þ trðW T WÞK1 W T U;
ð8Þ
123
852
Med Biol Eng Comput (2013) 51:849–858
images for recovering the geometry of the imaged object, an incandescent lamp was employed to provide epi-illumination. The exposure time for collecting fluorescent images with 2 9 2 binning was 2 s. 36 fluorescent images were obtained in 10° steps. 72 white light images were collected every 5° during a 360° rotation.
3 Results Several phantom experiments were employed to verify the proposed ATR method. In all the cases, the number of iterations was set as N = 8 and the stop regularization parameter was set as astop = 10-8 for the proposed ATR method. For comparisons, the conventional TR was also used in FMT reconstructions and the used regularization parameter is empirically chosen as a = 10-6, which can acquire the best reconstructed images in all the investigated regularization parameters (a = 10-3, 10-4, 10-5, 10-6, 10-7). Besides, the influence of the number of iterations N and the stop regularization parameter astop on the reconstruction were also investigated.
Fig. 2 Sections and 3D views of the reconstructed result for the single fluorescent target. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and reconstructed result using TR, respectively. All the images are normalized by the maximum of the results
3.1 Single fluorescent target without background fluorescence In this experiment, a transparent glass cylinder with a diameter of 2.4 cm, height of 4.8 cm filled with 1 % intralipid (la=0.02 cm-1 and l0 s=10 cm-1) was employed as the homogeneous phantom. No ICG was added in the intralipid solution so that there is no backgroud fluorescence in the phantom. A small transparent glass tube (0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM was used as the fluorescent target. Figure 2 shows the sections and 3D views of reconstructed results for the single fluorescent target without background fluorescence. Although both methods can correctly localize the embedded fluorescent target, reconstructed results with more accurate size of the fluorescent target and less noise are obtained with the proposed ATR method. To provide more details, the profiles along the dotted lines in Fig. 2 are shown in Fig. 3. The contrast-to-noise ratio (CNR) [2], relative error (RE) and full width at half maximum value (FWHM) are used to quantitatively evaluate the reconstructed results. The CNR is defined as lROI lROB CNR ¼ ; ð10Þ ðxROI rROI þ xROB rROB Þ1=2 where lROI and lROB are the mean values of the object and background regions in the reconstructed images, respectively; rROI and rROB are the corresponding standard
123
Fig. 3 Profiles taken through the sections shown in Fig. 2 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR, respectively
deviations. xROI and xROB are the ratio of the object and background regions with respect to the whole area, respectively. The contrast of the regions with fluorescent target to the background can be quantified by CNR, and greater CNR values imply a better image quality. RE is defined as RE ¼
kX Xtrue k ; kXtrue k
ð11Þ
where Xtrue is the true distributions of the fluorescent target and X is the reconstructed result. RE can comprehensively quantify reconstruction accuracy. The CNR, RE, and FWHM of the reconstructed results for the single fluorescent target without background fluorescence are listed in Table 1. The CNR of the reconstructed results obtained by ATR is greater than that of TR.
Med Biol Eng Comput (2013) 51:849–858
853
Table 1 Quantitative analysis of the reconstructed results for the single fluorescent target without background fluorescence Method
CNR
RE
FWHM (cm)
ATR
9.27
0.68
0.38
TR
6.09
1.44
0.59
The RE and FWHM of the reconstructed results obtained by ATR are smaller than those of TR. 3.2 Single fluorescent target with background fluorescence The voxels with no fluorescent probe can be seen as the background, which will be penalized with large regularization parameters. However, weak fluorescence in the background may mislead the adaptive choice of regularization parameters during the iteration of ATR. So a phantom experiment considering single fluorescent target with background fluorescence is carried out to assess the performance of the proposed ATR method. A small transparent glass tube (0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM was used as the fluorescent target. A transparent glass cylinder with a diameter of 2.4 cm, height of 4.8 cm filled with 1 % intralipid was employed as the homogeneous phantom. ICG was added in the intralipid solution to mimic the weak background fluorescence, and the ICG concentration in the intralipid solution is 0.65 nM which is 1/2,000 of the fluorescent target. Figure 4 shows the comparison of captured fluorescent images with and without background fluorescence. Although only a low concentration of ICG is presented in the background, the distributions and intensities of fluorescent signals have prominent changes. When the background fluorescence exists, the shape of reconstructed fluorescent target is distorted and a torus artifact is emerged, as shown in Fig. 5. The torus artifact of the reconstructed result using ATR is slighter than that of TR. The size of the reconstructed fluorescent target using ATR is more accurate than that of TR. The profiles along the dotted lines in Fig. 5 are shown in Fig. 6. Quantitative analysis is listed in Table 2. The CNR of the reconstructed results obtained by ATR is greater than that by TR. The RE and FWHM obtained by ATR are smaller than those by TR. The results demonstrate that the background fluorescence brings artifact into the reconstructed images, and ATR can suppress this artifact more effectively than the TR method. 3.3 Two fluorescent targets at the same height In this study, a transparent glass cylinder with 2.4 cm in diameter, 3.8 cm in height filled with 1 % intralipid was employed as the homogeneous phantom and two glass tubes
Fig. 4 The captured fluorescent images at different projections. The upper row is the case with no background fluorescence and the lower row is the case with background fluorescence
Fig. 5 Sections and 3D views of the reconstructed results. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and reconstructed result using TR, respectively. All the images are normalized by the maximum of the results
(0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM were used as the fluorescent targets. The two glass tubes were immersed in the phantom at the same height with an edge-to-edge distance of 0.5 cm [one was at (-0.2, 0.3, 1.6 cm), the other at (0.2, -0.3, 1.6 cm)]. Figure 7 shows the sections and 3D views of the reconstructed result for two fluorescent targets at the same height. The proposed ATR method can accurately separate the two fluorescent targets, and the sizes of the reconstructed fluorescent targets are generally consistent with the true sizes. But TR can barely separate the two fluorescent targets and the sizes of the reconstructed targets are larger than the true sizes. Furthermore, the reconstructed results using TR have more noise than that of the proposed method. To provide more details, the profiles along the dotted lines in Fig. 7 are shown in Fig. 8. In Fig. 7, the reconstructed results obtained by ATR have a higher spatial resolution than that by TR. Quantitative analysis of the reconstructed results for two fluorescent targets at the same height is listed in Table 3. The CNR of the reconstructed results obtained by ATR is
123
854
Med Biol Eng Comput (2013) 51:849–858
Fig. 6 Profiles taken through the sections shown in Fig. 5 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR, respectively
Table 2 Quantitative analysis of the reconstructed results for the single fluorescent target with background fluorescence Method
CNR
RE
FWHM (cm)
ATR
7.27
0.96
0.27
TR
3.38
2.47
0.52
Fig. 8 Profiles taken through the sections shown in Fig. 7 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR, respectively
Table 3 Quantitative analysis of the reconstructed results for two fluorescent targets at the same height Method
CNR
RE
FWHM (cm) Target 1
Target 2
ATR
8.68
0.22
0.41
0.39
TR
3.58
0.74
0.62
0.56
A transparent glass cylinder with 2.4 cm in diameter, 6.3 cm in height filled with 1 % intralipid was employed as the homogeneous phantom. Two glass tubes (0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM were used as the fluorescent targets. The two glass tubes were immersed in the phantom at different heights (one at the height of 5.2 cm, the other at 1.2 cm). Figure 9 is the reconstructed result for this case. ATR acquires both the correct location and size of the fluorescent targets. However, TR acquires only the correct locations of the fluorescent targets. Quantitative analysis of the results is listed in Table 4.
Fig. 7 Sections and 3D views of the reconstructed result for two fluorescent targets at the same height. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and TR, respectively. All the images are normalized by the maximum of the results
greater than that by TR. The RE and FWHM obtained by ATR are smaller than those by TR. 3.4 Two fluorescent targets at different heights To examine the whole body imaging performance, this study considers two fluorescent targets located at different heights with a relatively far distance between each other.
123
3.5 Two fluorescent targets with different ICG concentrations In this experiment, two glass tubes (0.3 cm in inner diameter) with different ICG concentrations were employed as the fluorescent targets. One was filled with 10-lL ICG with the concentration of 1.3 lM, and the other one was filled with 10-lL ICG with the concentration of 2.6 lM. As shown in Fig. 10, ATR can more accurately obtain the sizes of the fluorescent targets with less noise than TR. For both the methods, the size of reconstructed fluorescent target with high fluorescence concentration is a slightly larger than that with the low one. TR can more accurately
Med Biol Eng Comput (2013) 51:849–858
855
Fig. 9 3D views and sections of the reconstructed result for two fluorescent targets at different heights. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and TR, respectively. All the images are normalized by the maximum of the results
Table 4 Quantitative analysis of the reconstructed results for two fluorescent targets at different heights Method
CNR
RE
Height (cm) Target 1
Target 2
ATR
7.47
0.44
5.1
1.3
TR
4.67
0.90
5.3
1.1
Fig. 11 Profiles taken through the sections shown in Fig. 10 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR respectively Table 5 Quantitative analysis of the reconstructed results for two fluorescent targets with different ICG concentrations Method
Fig. 10 Sections and 3D views of the reconstructed result for two fluorescent targets with different ICG concentration. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and TR, respectively. All the images are normalized by the maximum of the results
obtain the relative fluorescence concentration of the fluorescent targets. The profiles (Fig. 11) can clearly illustrate this point. Quantitative analysis of the results is listed in Table 5. 3.6 Influence of the number of iterations N The proposed method is an iteration process, and the final step uses the prior information provided by all results of the previous iterations with different regularization parameters. It means that the final result is the comprehensive effect of
CNR
RE
FWHM (cm)
Normalized value
Target 1
Target 2
Target 1
Target 2
ATR
6.98
0.46
0.44
0.36
0.95
0.71
TR
4.10
0.95
0.68
0.27
0.96
0.56
all the regularization parameters used in the iteration process. The more regularization parameters are used in the iteration processes, the more prior information is provided into the final result. So a large number of the regularization parameters will improve the quality of the reconstructed result, but costs more computation time. Once the stop regularization parameter astop is fixed, the number of iterations N will decide the number of regularization parameters used in the iteration processes. So we examine the reconstructed results using different N when the stop regularization parameter is astop = 10-8. Figure 12 shows the sections of reconstructed results using different number of iterations N for the single and double fluorescent target case, respectively. With the increase of N, the reconstructed result has higher spatial resolution and less noise. When N C 10, the improvements of the reconstructed results are not evident. So N = 10 is suggested to make a compromise between reconstruction quality and computation time. Figure 13 is the quantitative analysis of the reconstructed results shown in Fig. 12. As N increases, the CNR increases and the RE decreases in the both cases. They demonstrate the improvements of the reconstructed results with the increase of N.
123
856
Med Biol Eng Comput (2013) 51:849–858
Fig. 13 CNR and RE of reconstructed results using different numbers of iterations N. a The single fluorescent target case. b The double fluorescent targets case Fig. 12 Sections of reconstructed results using different number of iterations N. a For single fluorescent target case. b For double fluorescent targets case. All the images are normalized by the maximum of the results
3.7 Influence of the stop regularization parameter astop Each iteration of the proposed method is a TR with a spatial varying regularization parameter. The stop regularization parameter astop will decide the range of all the regularization parameters, especially the smallest regularization parameter used in the final iteration. So we also evaluate the reconstructed results using different astop when the number of iterations is fixed to N = 10. Figure 14 shows the sections of reconstructed results using different astop for the two cases of the single and double fluorescent targets. With the decrease of astop, the reconstructed results have higher spatial resolution. When astop \ 10-9, the sizes of the reconstructed fluorescent targets are slightly smaller than that of the true ones. Unlike TR, a small regularization parameter can improve thespatial resolution without producing an additional noise. That is, the effects of spatial varying regularization parameters based on the prior information provided by all results of the previous iterations. As a result, the proposed method can use a smaller regularization parameter than the TR to improve the spatial resolution, without introducing additional noise. Figure 15 is the quantitative analysis of the reconstructed results shown in Fig. 14. When astop C 10-9, with the decrease of astop, the CNR increases and the RE decreases in both cases. The improvements of the
123
Fig. 14 Sections of reconstructed results using different stop regularization parameters astop. a The single fluorescent target case. b The double fluorescent targets case. All the images are normalized by the maximum of the results
reconstructed results come from the improvements of spatial resolution due to the decrease of astop. When astop B 10-9, the CNR decreases and the RE increases in both cases. These results demonstrate that when the regularization parameters are too small, the reconstructed
Med Biol Eng Comput (2013) 51:849–858
Fig. 15 CNR and RE of reconstructed results using different stop regularization parameter astop. a The single fluorescent target case. b The double fluorescent targets case
results will be degenerated. So the optimal range of astop is 10-9 B astop B 10-8.
4 Discussion Due to the highly ill-posed nature of the inverse problem in FMT, the reconstruction is vulnerable to noise in the measurements. Regularization technique is a powerful tool to obtain biologically reasonable solutions. Regularization parameter is a key factor which determines the spatial resolution and noise of the reconstructed result. A too large parameter makes the reconstructed image over-smooth with low spatial resolution and less noise, while a too small parameter makes the reconstructed image over-sharp with high spatial resolution and more noise. Generally, the optimal regularization parameter makes a compromise between the spatial resolution and noise of the reconstructed image, which usually leads to a poor spatial resolution of the reconstructed image for noise suppression. The proposed ATR method uses different regularization parameters to spatially penalize different voxels of the reconstructed image. The larger regularization parameters are used to spatially penalize the voxels which have no fluorescent target distributions and the smaller regularization parameters are used to penalize the voxels which have fluorescent target distributions. The spatially varying regularization parameters are extracted from the results of the previous iterations without complementary information
857
provided by other imaging modalities. The regularization parameters for region with fluorescent target distributions are smaller than the optimal regularization parameter used in TR, and the regularization parameters for region with fluorescent target distributions are larger than the optimal regularization parameter used in TR. So the proposed ATR method can adequately sharpen or smoothen the region with or without fluorescent target distributions, respectively. As a result, the reconstructed results by the proposed method have a higher spatial resolution and less noises compared with TR. In the proposed ATR method, the number of iterations N and the stop regularization parameter astop jointly decide the spatially varying regularization parameters in all the iterations. With the increase of N, the reconstructed results are improved at the cost of computation time and the improvements of the reconstructed results are not significant when N C 10. So N = 10 is suggested to avoid too heavy computation burden. Too large or too small value of astop will lead to over-smooth or over-sharpen reconstructed results, respectively. So the suggested range of values for astop is 10-9,10-8. These suggested N and astop are probably not optimal for data obtained by different experimental setups in different application circumstances, but we hope they can provide significative references. The results of phantom experiments demonstrate that the proposed ATR method can improve the spatial resolution and reduce the noise of FMT reconstruction at the same time. The ATR method is an iteration algorithm and each iteration is a TR method with spatially varying regularization parameters. So the proposed ATR method costs more time than TR. However, this speed is acceptable without considering real-time reconstruction in most situations. 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; the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114; National Science and technology support program under Grant No. 2012BAI23B00 and 2011BAI02B03.
References 1. Axelsson J, Svensson S (2007) Spatially varying regularization based on spectrally resolved fluorescence emission in fluorescence molecular tomography. Opt Express 15(21):13574–13584 2. Baritaux JC, Hassler, Bucher KM, Sanyal S, Unser M (2011) Sparsity-driven reconstruction for FDOT with anatomical priors. IEEE Trans Med Imaging 30(5):1143–1153 3. Ducros N, D’andrea C, Valentini G, Rudge T, Arridge S (2010) Bassi A Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination. Opt Lett 35(21):3676– 3678
123
858 4. Ducros N, Bassi A, Valentini G, Schweiger M, Arridge S, D’Andrea C (2011) Multiple-view fluorescence optical tomography reconstruction using compression of experimental data. Opt Lett 36(8):1377–1379 5. Golub GH, Heath M, Wahba G (1979) Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21(2):215–223 6. Graves EE, Weissleder R, Ntziachristos V (2004) Fluorescence molecular imaging of small animal tumor models. Curr Mol Med 4(4): 419–430 7. Hansen PC, O’Leary DP (1993) The use of the l-curve in the regularization of discrete ill-posed problems. SIAM J Sci Comput 14(6):1487–1503 8. Hansen PC (1992) Analysis of discrete ill-posed problems by means of the l-curve. SIAM Rev 34(4):561–580 9. Hiltunen P, Calvetti D, Somersalo E (2008) An adaptive smoothness regularization algorithm for optical tomography. Opt Express 16(24):19957–19977 10. Hyde D, Miller E, Brooks D, Ntziachristos V (2010) Data specific spatially varying regularization for multimodal fluorescence molecular tomography. IEEE Trans Med Imaging 29(2):365–374 11. Katamreddy SH, Yalavarthy PK (2012) Model-resolution based regularization improves near infrared diffuse optical tomography. J Opt Soc Am A 29(5):649–656 12. Liu F, Liu X, Wang DF, Zhang B, Bai J (2010) A parallel excitation based fluorescence molecular tomography system for whole-body simultaneous imaging of small animals. Ann Biomed Eng 38(11):3440–3448 13. Massoud TF, Gambhir SS (2003) Molecular imaging in living subjects: seeing fundamental biological processes in a new light. Gene Dev 17(5):545–580 14. Montet X, Ntziachristos V, Grimm J, Weissleder R (2005) Tomographic fluorescence mapping of tumor targets. Cancer Res 65(14):6330–6336 15. Niu H, Guo P, Ji L, Zhao Q, Jiang T (2008) Improving image quality of diffuse optical tomography with a projection-error-
123
Med Biol Eng Comput (2013) 51:849–858
16.
17.
18.
19.
20.
21. 22.
23.
24.
25.
26.
based adaptive regularization method. Opt Express 16(17): 12423–12434 Ntziachristos V, Bremer C, Graves EE, Ripoll J, Weissleder R (2002) In vivo tomographic imaging of near-infrared fluorescent probes. Mol imaging 1(2):82–88 Pogue BW, McBride TO, Prewitt J, O’sterberg UL, Paulsen KD (1999) Spatially variant regularization improves diffuse optical tomography. Appl Opt 38(13):2950–2961 Ripoll J, Ntziachristos V, Carminati R, Nieto-Vesperinas M (2001) Kirchhoff approximation for diffusive waves. Phys Rev E 64:051917 Ripoll J, Nieto-Vesperinas M, Weissleder R, Ntziachristos V (2002) Fast analytical approximation for arbitrary geometries in diffuse optical tomography. Opt Lett 27(7):527–529 Rudge TJ, Soloviev VY, Arridge SR (2010) Fast image reconstruction in fluoresence optical tomography using data compression. Opt Lett 35(5):763–765 Rudin M, Weissleder R (2003) Molecular imaging in drug discovery and development. Nat Rev Drug Discov 2(2):123–131 Scherzer O (1993) The use of morozov’s discrepancy principle for tikhonov regularization for solving nonlinear ill-posed problems. Comput Lett 51:45–60 Schulz RB, Ripoll J, Ntziachristos V (2004) Experimental fluorescence tomography of tissues with noncontact measurements. IEEE Trans Med Imaging 23(4):492–500 Schweiger M, Arridge SR, Hiraoka M, Delpy DT (1995) The finite-element method for the propagation of light in scattering media—boundary and source conditions. Med Phys 22(11):1779– 1792 Srinivasan S, Pogue B, Dehghani H, Jiang S, Song X, Paulsen K (2004) Improved quantification of small objects in near-infrared diffuse optical tomography. J Biomed Opt 9(6):1161–1171 Tikhonov AN, Arsenin VY (1977) Solutions of ill-posed problems. Wiley, Washington, DC
ORIGINAL ARTICLE
An adaptive Tikhonov regularization method for fluorescence molecular tomography Xu Cao • Bin Zhang • Xin Wang • Fei Liu Ke Liu • Jianwen Luo • Jing Bai
•
Received: 8 November 2012 / Accepted: 23 February 2013 / Published online: 16 March 2013 Ó International Federation for Medical and Biological Engineering 2013
Abstract The high degree of absorption and scattering of photons propagating through biological tissues makes fluorescence molecular tomography (FMT) reconstruction a severe ill-posed problem and the reconstructed result is susceptible to noise in the measurements. To obtain a reasonable solution, Tikhonov regularization (TR) is generally employed to solve the inverse problem of FMT. However, with a fixed regularization parameter, the Tikhonov solutions suffer from low resolution. In this work, an adaptive Tikhonov regularization (ATR) method is presented. Considering that large regularization parameters can smoothen the solution with low spatial resolution, while small regularization parameters can sharpen the solution with high level of noise, the ATR method adaptively updates the spatially varying regularization parameters during the iteration process and uses them to penalize the solutions. The ATR method can adequately sharpen the feasible region with fluorescent probes and smoothen the region without fluorescent probes resorting to no complementary priori information. Phantom experiments are performed to verify the feasibility of the proposed method. The results demonstrate that the proposed method can X. Cao B. Zhang X. Wang F. Liu J. Bai (&) Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China e-mail: [email protected]; [email protected] K. Liu Division of Electronics and Information Technology, National Institute of Metrology, Beijing 100013, People’s Republic of China J. Luo Center for Biomedical Imaging Research and Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China
improve the spatial resolution and reduce the noise of FMT reconstruction at the same time. Keywords Adaptive regularization Fluorescence molecular tomography regularization Reconstruction Tikhonov
1 Introduction Fluorescence molecular tomography (FMT) can reveal biological processes in vivo by three-dimensional (3D) visualization of specific fluorescent probes [16]. With the advantages of low cost, non-invasive and non-ionizing, FMT has been developed into a powerful in vivo preclinical imaging method for studies of cancer diagnosis [6, 14], drug discovery [21], and gene expression visualization [13]. As a typical inverse problem, the high degree of absorption and scattering of photons propagating through biological tissues makes the FMT reconstruction a severe ill-posed problem. It means that the solution is susceptible to noise in the measurements. Tikhonov regularization (TR) is a widespread method to deal with the ill-posed inverse problems, in which an L2-norm constraint term is added into the data-fitting term to improve the stability [26]. Regularization parameter is important for regularization method in controlling the smooth degree of the reconstructed image. Generally, a large regularization parameter results in a smooth image with low spatial resolution, while a small regularization parameter results in high spatial resolution image with lots of high-frequency noise. To make a compromise between the spatial resolution and signal-to-noise ratio of the reconstructed image, several methods are developed to choose an optimal
123
850
regularization parameter, such as the L-curve method [8, 7], the cross-validation method [5], and the discrepancy principle [22]. However, in practice, the selected regularization parameter usually leads to a poor spatial resolution in the reconstructed FMT images, and to make matters worse, the selected regularization parameter is not the optimal one due to the serious ill-posed nature and the model mismatch of FMT. So the optimal regularization parameter is often empirically chosen in practical applications [4, 3, 20]. To improve the quality of the reconstructed image using single constant regularization parameter, several adaptive regularization methods are proposed to improve the resolution and contrast of the reconstructed image in diffuse optical tomography (DOT) [9, 11, 15, 17, 25]. Similarly, the use of spatially varying regularization parameters according to the additional information obtained from multi-spectral fluorescent emission recordings [1] or the prior anatomical information obtained from CT images [10] in the reconstruction has been reported to improve image quality in FMT. But these two methods need complementary information which increases the complexity of imaging system. Due to the differences between FMT and DOT, the adaptive regularization methods proposed for DOT reconstruction cannot be applied in FMT reconstruction, and to the best of our knowledge, there is not an appropriate adaptive regularization method independent of complementary information for FMT reconstruction. In this study, we modify the conventional TR and propose an adaptive Tikhonov regularization (ATR) without using the complementary information, which dynamically updates the spatially varying regularization parameters based on results of previous iterations. The principle of this method is using different regularization parameters to spatially penalize different voxels of the reconstructed image. Large regularization parameters are used to spatially penalize the voxels which have no fluorescent probes and small regularization parameters are used to penalize the voxels which have fluorescent probes. As a result, the region with fluorescent probes can be adequately sharpened and the region without fluorescent probes can be adequately smoothened. To adaptively update the regularization parameters, the iteration starts with a large regularization parameter which is much larger than the empirically chosen one. During iteration processes, result of each iteration is used as the prior information to decide the changes of the regularization parameters in the next iteration. The regularization parameters corresponding to the zero value voxels stay the same and those corresponding to the non-zero value voxels are reduced. So the region with fluorescent probes is sharpened by smaller and smaller regularization parameters and the region without fluorescent target distribution is smoothened by large regularization parameters. Finally a
123
Med Biol Eng Comput (2013) 51:849–858
robust FMT reconstruction with high spatial resolution can be acquired. To examine the improvements of spatial resolution and reduction of the noise in FMT reconstruction, phantom experiments are implemented, respectively. Furthermore, the influences of important parameters of the proposed method on the reconstructed image are discussed.
2 Methods 2.1 Forward model Considering a continuous wave (CW) point excitation source at rs, the diffusion equation with the Robin-type boundary condition [24] to describe the light transportation field in tissue is r DrGðrs ; rÞ la Gðrs ; rÞ ¼ dðr rs Þ;
ð1Þ
where r is the position vector in the domain, b is the position vector on the boundary, la is the absorption coefficient, D = 1/(3(l a ? l0 s)) is the diffusion coefficient of the tissue with the reduced scattering coefficient l0 s. 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/l0 s into the medium from the illumination spot. The Robin boundary is given by 2DðbÞ
oGðrs ; bÞ þ qGðrs ; bÞ ¼ 0; on
ð2Þ
where 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. The fluorescent measurements /(rs, rd) detected at location rd due to an illumination spot located at rs can be formulated as follows, using the first Rytov approximation [23] Z /ðrs ; rd Þ ¼ H Gðrs ; rÞxðrÞGðrd ; rÞd3 r; ð3Þ where G(rs, r) and G(rd, r) are the Green’s functions of excitation and emission light respectively, H is a unit-less constant taking account of the unknown gain and attenuation factors of the system. x(r) is the distribution of fluorescent probe representing the fluorochrome concentration at location r multiplied by the fluorescent yield. After the imaged object is discretized, a Kirchhoff Approximation (KA) [19, 18] solution of Eq. (1) is substituted into Eq. (3) /ðrs ; rd Þ ¼ ws;d X;
ð4Þ
where the column vector X is the discrete formation of x(r) representing the distribution of fluorescent probe to be
Med Biol Eng Comput (2013) 51:849–858
851
reconstructed, ws,d is the source detector map with the element of ws;d ðrÞ ¼ HDvGðrs ; rÞGðrd ; rÞ;
pffi computed by c ¼ ½Nastop and N is the number of iterations. The reconstruction strategies can be summarized as follows
ð5Þ
where Dv is the volume of the discrete voxel. Finally, the following matrix equation of the inverse problem is formed by assembling all the source detector maps as follows U ¼ WX;
ð6Þ
where W = {ws,d} is the weight matrix mapping unknown distributions of the fluorescent probes inside the small animal onto the measured fluorescent data over the surface. U ¼ f/ðrs ; rd Þg is a column vector constituted by fluorescent measurements.
2.3 Experimental setup
where K is a diagonal matrix with diagonal elements of an(r) which is updated according to the previous iteration result Xn-1. The negative elements of Xn-1 are physically impossible, so those are reset to zero. an1 ðrÞ Xn1 ðrÞ ¼ 0 an ðrÞ ¼ ð9Þ can1 ðrÞ Xn1 ðrÞ [ 0
Phantom experiments were conducted based on the FMT system previously developed in our laboratory [12]. The schematic diagram of the system is shown in Fig. 1. The imaged object was fixed onto the rotation stage, which can rotate around its longitudinal axis for full-angle projection acquisitions. An appropriate rotation speed (6°/s) was selected to avoid skew and internal organs movement during the in vivo imaging process. The point excitation source was provided by a 250 W Halogen lamp (7ILT250, 7-star, Beijing, China). The excitation light traveled through a 775 ± 23 nm 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 ± 6 nm 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
where 0 \ c \ 1 is a shrinkage factor of the regularization parameter. In this equation, an(r) is spatially varying and adaptively updated based on the previous iteration result Xn-1. an(r) corresponding to zero value elements (i.e., voxels have no fluorescent probes) of Xn-1 stay the same, while an(r) corresponding to the non-zero value ones (i.e., voxels have fluorescent probes) are reduced. In the next iteration, a larger regularization parameter is used to spatially penalize the voxels which have no probes and a small one is used to spatially penalize the voxels which have probes. The initial regularization parameter is set as a0 = 1. So the first solution is equivalent to Eq. (7). During the iterations, the regularization parameters corresponding to the voxels with fluorescent probes get smaller and smaller. To avoid unstable solutions caused by too small regularization parameters, stop regularization parameter is set as astop. So the shrinkage factor can be
Fig. 1 The schematic of the experimental setup
2.2 Adaptive Tikhonov regularization method The matrix W on the right of Eq. (6) is generally ill-posed and TR is used to obtain a stable solution X ¼ ½W T W þ trðW T WÞaI1 W T U;
ð7Þ
where I is the identity matrix, tr(*) is trace of a matrix and a is the regularization parameter to stabilize the solution. In this study, Eq. (7) is modified into an iteration process using an adaptive spatially varying regularization parameter an(r) rather than a single scalar value a Xn ¼ ½W T W þ trðW T WÞK1 W T U;
ð8Þ
123
852
Med Biol Eng Comput (2013) 51:849–858
images for recovering the geometry of the imaged object, an incandescent lamp was employed to provide epi-illumination. The exposure time for collecting fluorescent images with 2 9 2 binning was 2 s. 36 fluorescent images were obtained in 10° steps. 72 white light images were collected every 5° during a 360° rotation.
3 Results Several phantom experiments were employed to verify the proposed ATR method. In all the cases, the number of iterations was set as N = 8 and the stop regularization parameter was set as astop = 10-8 for the proposed ATR method. For comparisons, the conventional TR was also used in FMT reconstructions and the used regularization parameter is empirically chosen as a = 10-6, which can acquire the best reconstructed images in all the investigated regularization parameters (a = 10-3, 10-4, 10-5, 10-6, 10-7). Besides, the influence of the number of iterations N and the stop regularization parameter astop on the reconstruction were also investigated.
Fig. 2 Sections and 3D views of the reconstructed result for the single fluorescent target. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and reconstructed result using TR, respectively. All the images are normalized by the maximum of the results
3.1 Single fluorescent target without background fluorescence In this experiment, a transparent glass cylinder with a diameter of 2.4 cm, height of 4.8 cm filled with 1 % intralipid (la=0.02 cm-1 and l0 s=10 cm-1) was employed as the homogeneous phantom. No ICG was added in the intralipid solution so that there is no backgroud fluorescence in the phantom. A small transparent glass tube (0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM was used as the fluorescent target. Figure 2 shows the sections and 3D views of reconstructed results for the single fluorescent target without background fluorescence. Although both methods can correctly localize the embedded fluorescent target, reconstructed results with more accurate size of the fluorescent target and less noise are obtained with the proposed ATR method. To provide more details, the profiles along the dotted lines in Fig. 2 are shown in Fig. 3. The contrast-to-noise ratio (CNR) [2], relative error (RE) and full width at half maximum value (FWHM) are used to quantitatively evaluate the reconstructed results. The CNR is defined as lROI lROB CNR ¼ ; ð10Þ ðxROI rROI þ xROB rROB Þ1=2 where lROI and lROB are the mean values of the object and background regions in the reconstructed images, respectively; rROI and rROB are the corresponding standard
123
Fig. 3 Profiles taken through the sections shown in Fig. 2 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR, respectively
deviations. xROI and xROB are the ratio of the object and background regions with respect to the whole area, respectively. The contrast of the regions with fluorescent target to the background can be quantified by CNR, and greater CNR values imply a better image quality. RE is defined as RE ¼
kX Xtrue k ; kXtrue k
ð11Þ
where Xtrue is the true distributions of the fluorescent target and X is the reconstructed result. RE can comprehensively quantify reconstruction accuracy. The CNR, RE, and FWHM of the reconstructed results for the single fluorescent target without background fluorescence are listed in Table 1. The CNR of the reconstructed results obtained by ATR is greater than that of TR.
Med Biol Eng Comput (2013) 51:849–858
853
Table 1 Quantitative analysis of the reconstructed results for the single fluorescent target without background fluorescence Method
CNR
RE
FWHM (cm)
ATR
9.27
0.68
0.38
TR
6.09
1.44
0.59
The RE and FWHM of the reconstructed results obtained by ATR are smaller than those of TR. 3.2 Single fluorescent target with background fluorescence The voxels with no fluorescent probe can be seen as the background, which will be penalized with large regularization parameters. However, weak fluorescence in the background may mislead the adaptive choice of regularization parameters during the iteration of ATR. So a phantom experiment considering single fluorescent target with background fluorescence is carried out to assess the performance of the proposed ATR method. A small transparent glass tube (0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM was used as the fluorescent target. A transparent glass cylinder with a diameter of 2.4 cm, height of 4.8 cm filled with 1 % intralipid was employed as the homogeneous phantom. ICG was added in the intralipid solution to mimic the weak background fluorescence, and the ICG concentration in the intralipid solution is 0.65 nM which is 1/2,000 of the fluorescent target. Figure 4 shows the comparison of captured fluorescent images with and without background fluorescence. Although only a low concentration of ICG is presented in the background, the distributions and intensities of fluorescent signals have prominent changes. When the background fluorescence exists, the shape of reconstructed fluorescent target is distorted and a torus artifact is emerged, as shown in Fig. 5. The torus artifact of the reconstructed result using ATR is slighter than that of TR. The size of the reconstructed fluorescent target using ATR is more accurate than that of TR. The profiles along the dotted lines in Fig. 5 are shown in Fig. 6. Quantitative analysis is listed in Table 2. The CNR of the reconstructed results obtained by ATR is greater than that by TR. The RE and FWHM obtained by ATR are smaller than those by TR. The results demonstrate that the background fluorescence brings artifact into the reconstructed images, and ATR can suppress this artifact more effectively than the TR method. 3.3 Two fluorescent targets at the same height In this study, a transparent glass cylinder with 2.4 cm in diameter, 3.8 cm in height filled with 1 % intralipid was employed as the homogeneous phantom and two glass tubes
Fig. 4 The captured fluorescent images at different projections. The upper row is the case with no background fluorescence and the lower row is the case with background fluorescence
Fig. 5 Sections and 3D views of the reconstructed results. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and reconstructed result using TR, respectively. All the images are normalized by the maximum of the results
(0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM were used as the fluorescent targets. The two glass tubes were immersed in the phantom at the same height with an edge-to-edge distance of 0.5 cm [one was at (-0.2, 0.3, 1.6 cm), the other at (0.2, -0.3, 1.6 cm)]. Figure 7 shows the sections and 3D views of the reconstructed result for two fluorescent targets at the same height. The proposed ATR method can accurately separate the two fluorescent targets, and the sizes of the reconstructed fluorescent targets are generally consistent with the true sizes. But TR can barely separate the two fluorescent targets and the sizes of the reconstructed targets are larger than the true sizes. Furthermore, the reconstructed results using TR have more noise than that of the proposed method. To provide more details, the profiles along the dotted lines in Fig. 7 are shown in Fig. 8. In Fig. 7, the reconstructed results obtained by ATR have a higher spatial resolution than that by TR. Quantitative analysis of the reconstructed results for two fluorescent targets at the same height is listed in Table 3. The CNR of the reconstructed results obtained by ATR is
123
854
Med Biol Eng Comput (2013) 51:849–858
Fig. 6 Profiles taken through the sections shown in Fig. 5 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR, respectively
Table 2 Quantitative analysis of the reconstructed results for the single fluorescent target with background fluorescence Method
CNR
RE
FWHM (cm)
ATR
7.27
0.96
0.27
TR
3.38
2.47
0.52
Fig. 8 Profiles taken through the sections shown in Fig. 7 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR, respectively
Table 3 Quantitative analysis of the reconstructed results for two fluorescent targets at the same height Method
CNR
RE
FWHM (cm) Target 1
Target 2
ATR
8.68
0.22
0.41
0.39
TR
3.58
0.74
0.62
0.56
A transparent glass cylinder with 2.4 cm in diameter, 6.3 cm in height filled with 1 % intralipid was employed as the homogeneous phantom. Two glass tubes (0.3 cm in inner diameter) filled with 10-lL ICG with the concentration of 1.3 lM were used as the fluorescent targets. The two glass tubes were immersed in the phantom at different heights (one at the height of 5.2 cm, the other at 1.2 cm). Figure 9 is the reconstructed result for this case. ATR acquires both the correct location and size of the fluorescent targets. However, TR acquires only the correct locations of the fluorescent targets. Quantitative analysis of the results is listed in Table 4.
Fig. 7 Sections and 3D views of the reconstructed result for two fluorescent targets at the same height. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and TR, respectively. All the images are normalized by the maximum of the results
greater than that by TR. The RE and FWHM obtained by ATR are smaller than those by TR. 3.4 Two fluorescent targets at different heights To examine the whole body imaging performance, this study considers two fluorescent targets located at different heights with a relatively far distance between each other.
123
3.5 Two fluorescent targets with different ICG concentrations In this experiment, two glass tubes (0.3 cm in inner diameter) with different ICG concentrations were employed as the fluorescent targets. One was filled with 10-lL ICG with the concentration of 1.3 lM, and the other one was filled with 10-lL ICG with the concentration of 2.6 lM. As shown in Fig. 10, ATR can more accurately obtain the sizes of the fluorescent targets with less noise than TR. For both the methods, the size of reconstructed fluorescent target with high fluorescence concentration is a slightly larger than that with the low one. TR can more accurately
Med Biol Eng Comput (2013) 51:849–858
855
Fig. 9 3D views and sections of the reconstructed result for two fluorescent targets at different heights. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and TR, respectively. All the images are normalized by the maximum of the results
Table 4 Quantitative analysis of the reconstructed results for two fluorescent targets at different heights Method
CNR
RE
Height (cm) Target 1
Target 2
ATR
7.47
0.44
5.1
1.3
TR
4.67
0.90
5.3
1.1
Fig. 11 Profiles taken through the sections shown in Fig. 10 indicated by the dotted lines. The solid line, the solid line with star markers and the solid line with circle markers are the profiles of the true fluorescent target, reconstructed result using ATR and TR respectively Table 5 Quantitative analysis of the reconstructed results for two fluorescent targets with different ICG concentrations Method
Fig. 10 Sections and 3D views of the reconstructed result for two fluorescent targets with different ICG concentration. The first, second and third columns are the true fluorescent target, reconstructed result using ATR and TR, respectively. All the images are normalized by the maximum of the results
obtain the relative fluorescence concentration of the fluorescent targets. The profiles (Fig. 11) can clearly illustrate this point. Quantitative analysis of the results is listed in Table 5. 3.6 Influence of the number of iterations N The proposed method is an iteration process, and the final step uses the prior information provided by all results of the previous iterations with different regularization parameters. It means that the final result is the comprehensive effect of
CNR
RE
FWHM (cm)
Normalized value
Target 1
Target 2
Target 1
Target 2
ATR
6.98
0.46
0.44
0.36
0.95
0.71
TR
4.10
0.95
0.68
0.27
0.96
0.56
all the regularization parameters used in the iteration process. The more regularization parameters are used in the iteration processes, the more prior information is provided into the final result. So a large number of the regularization parameters will improve the quality of the reconstructed result, but costs more computation time. Once the stop regularization parameter astop is fixed, the number of iterations N will decide the number of regularization parameters used in the iteration processes. So we examine the reconstructed results using different N when the stop regularization parameter is astop = 10-8. Figure 12 shows the sections of reconstructed results using different number of iterations N for the single and double fluorescent target case, respectively. With the increase of N, the reconstructed result has higher spatial resolution and less noise. When N C 10, the improvements of the reconstructed results are not evident. So N = 10 is suggested to make a compromise between reconstruction quality and computation time. Figure 13 is the quantitative analysis of the reconstructed results shown in Fig. 12. As N increases, the CNR increases and the RE decreases in the both cases. They demonstrate the improvements of the reconstructed results with the increase of N.
123
856
Med Biol Eng Comput (2013) 51:849–858
Fig. 13 CNR and RE of reconstructed results using different numbers of iterations N. a The single fluorescent target case. b The double fluorescent targets case Fig. 12 Sections of reconstructed results using different number of iterations N. a For single fluorescent target case. b For double fluorescent targets case. All the images are normalized by the maximum of the results
3.7 Influence of the stop regularization parameter astop Each iteration of the proposed method is a TR with a spatial varying regularization parameter. The stop regularization parameter astop will decide the range of all the regularization parameters, especially the smallest regularization parameter used in the final iteration. So we also evaluate the reconstructed results using different astop when the number of iterations is fixed to N = 10. Figure 14 shows the sections of reconstructed results using different astop for the two cases of the single and double fluorescent targets. With the decrease of astop, the reconstructed results have higher spatial resolution. When astop \ 10-9, the sizes of the reconstructed fluorescent targets are slightly smaller than that of the true ones. Unlike TR, a small regularization parameter can improve thespatial resolution without producing an additional noise. That is, the effects of spatial varying regularization parameters based on the prior information provided by all results of the previous iterations. As a result, the proposed method can use a smaller regularization parameter than the TR to improve the spatial resolution, without introducing additional noise. Figure 15 is the quantitative analysis of the reconstructed results shown in Fig. 14. When astop C 10-9, with the decrease of astop, the CNR increases and the RE decreases in both cases. The improvements of the
123
Fig. 14 Sections of reconstructed results using different stop regularization parameters astop. a The single fluorescent target case. b The double fluorescent targets case. All the images are normalized by the maximum of the results
reconstructed results come from the improvements of spatial resolution due to the decrease of astop. When astop B 10-9, the CNR decreases and the RE increases in both cases. These results demonstrate that when the regularization parameters are too small, the reconstructed
Med Biol Eng Comput (2013) 51:849–858
Fig. 15 CNR and RE of reconstructed results using different stop regularization parameter astop. a The single fluorescent target case. b The double fluorescent targets case
results will be degenerated. So the optimal range of astop is 10-9 B astop B 10-8.
4 Discussion Due to the highly ill-posed nature of the inverse problem in FMT, the reconstruction is vulnerable to noise in the measurements. Regularization technique is a powerful tool to obtain biologically reasonable solutions. Regularization parameter is a key factor which determines the spatial resolution and noise of the reconstructed result. A too large parameter makes the reconstructed image over-smooth with low spatial resolution and less noise, while a too small parameter makes the reconstructed image over-sharp with high spatial resolution and more noise. Generally, the optimal regularization parameter makes a compromise between the spatial resolution and noise of the reconstructed image, which usually leads to a poor spatial resolution of the reconstructed image for noise suppression. The proposed ATR method uses different regularization parameters to spatially penalize different voxels of the reconstructed image. The larger regularization parameters are used to spatially penalize the voxels which have no fluorescent target distributions and the smaller regularization parameters are used to penalize the voxels which have fluorescent target distributions. The spatially varying regularization parameters are extracted from the results of the previous iterations without complementary information
857
provided by other imaging modalities. The regularization parameters for region with fluorescent target distributions are smaller than the optimal regularization parameter used in TR, and the regularization parameters for region with fluorescent target distributions are larger than the optimal regularization parameter used in TR. So the proposed ATR method can adequately sharpen or smoothen the region with or without fluorescent target distributions, respectively. As a result, the reconstructed results by the proposed method have a higher spatial resolution and less noises compared with TR. In the proposed ATR method, the number of iterations N and the stop regularization parameter astop jointly decide the spatially varying regularization parameters in all the iterations. With the increase of N, the reconstructed results are improved at the cost of computation time and the improvements of the reconstructed results are not significant when N C 10. So N = 10 is suggested to avoid too heavy computation burden. Too large or too small value of astop will lead to over-smooth or over-sharpen reconstructed results, respectively. So the suggested range of values for astop is 10-9,10-8. These suggested N and astop are probably not optimal for data obtained by different experimental setups in different application circumstances, but we hope they can provide significative references. The results of phantom experiments demonstrate that the proposed ATR method can improve the spatial resolution and reduce the noise of FMT reconstruction at the same time. The ATR method is an iteration algorithm and each iteration is a TR method with spatially varying regularization parameters. So the proposed ATR method costs more time than TR. However, this speed is acceptable without considering real-time reconstruction in most situations. 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; the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114; National Science and technology support program under Grant No. 2012BAI23B00 and 2011BAI02B03.
References 1. Axelsson J, Svensson S (2007) Spatially varying regularization based on spectrally resolved fluorescence emission in fluorescence molecular tomography. Opt Express 15(21):13574–13584 2. Baritaux JC, Hassler, Bucher KM, Sanyal S, Unser M (2011) Sparsity-driven reconstruction for FDOT with anatomical priors. IEEE Trans Med Imaging 30(5):1143–1153 3. Ducros N, D’andrea C, Valentini G, Rudge T, Arridge S (2010) Bassi A Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination. Opt Lett 35(21):3676– 3678
123
858 4. Ducros N, Bassi A, Valentini G, Schweiger M, Arridge S, D’Andrea C (2011) Multiple-view fluorescence optical tomography reconstruction using compression of experimental data. Opt Lett 36(8):1377–1379 5. Golub GH, Heath M, Wahba G (1979) Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21(2):215–223 6. Graves EE, Weissleder R, Ntziachristos V (2004) Fluorescence molecular imaging of small animal tumor models. Curr Mol Med 4(4): 419–430 7. Hansen PC, O’Leary DP (1993) The use of the l-curve in the regularization of discrete ill-posed problems. SIAM J Sci Comput 14(6):1487–1503 8. Hansen PC (1992) Analysis of discrete ill-posed problems by means of the l-curve. SIAM Rev 34(4):561–580 9. Hiltunen P, Calvetti D, Somersalo E (2008) An adaptive smoothness regularization algorithm for optical tomography. Opt Express 16(24):19957–19977 10. Hyde D, Miller E, Brooks D, Ntziachristos V (2010) Data specific spatially varying regularization for multimodal fluorescence molecular tomography. IEEE Trans Med Imaging 29(2):365–374 11. Katamreddy SH, Yalavarthy PK (2012) Model-resolution based regularization improves near infrared diffuse optical tomography. J Opt Soc Am A 29(5):649–656 12. Liu F, Liu X, Wang DF, Zhang B, Bai J (2010) A parallel excitation based fluorescence molecular tomography system for whole-body simultaneous imaging of small animals. Ann Biomed Eng 38(11):3440–3448 13. Massoud TF, Gambhir SS (2003) Molecular imaging in living subjects: seeing fundamental biological processes in a new light. Gene Dev 17(5):545–580 14. Montet X, Ntziachristos V, Grimm J, Weissleder R (2005) Tomographic fluorescence mapping of tumor targets. Cancer Res 65(14):6330–6336 15. Niu H, Guo P, Ji L, Zhao Q, Jiang T (2008) Improving image quality of diffuse optical tomography with a projection-error-
123
Med Biol Eng Comput (2013) 51:849–858
16.
17.
18.
19.
20.
21. 22.
23.
24.
25.
26.
based adaptive regularization method. Opt Express 16(17): 12423–12434 Ntziachristos V, Bremer C, Graves EE, Ripoll J, Weissleder R (2002) In vivo tomographic imaging of near-infrared fluorescent probes. Mol imaging 1(2):82–88 Pogue BW, McBride TO, Prewitt J, O’sterberg UL, Paulsen KD (1999) Spatially variant regularization improves diffuse optical tomography. Appl Opt 38(13):2950–2961 Ripoll J, Ntziachristos V, Carminati R, Nieto-Vesperinas M (2001) Kirchhoff approximation for diffusive waves. Phys Rev E 64:051917 Ripoll J, Nieto-Vesperinas M, Weissleder R, Ntziachristos V (2002) Fast analytical approximation for arbitrary geometries in diffuse optical tomography. Opt Lett 27(7):527–529 Rudge TJ, Soloviev VY, Arridge SR (2010) Fast image reconstruction in fluoresence optical tomography using data compression. Opt Lett 35(5):763–765 Rudin M, Weissleder R (2003) Molecular imaging in drug discovery and development. Nat Rev Drug Discov 2(2):123–131 Scherzer O (1993) The use of morozov’s discrepancy principle for tikhonov regularization for solving nonlinear ill-posed problems. Comput Lett 51:45–60 Schulz RB, Ripoll J, Ntziachristos V (2004) Experimental fluorescence tomography of tissues with noncontact measurements. IEEE Trans Med Imaging 23(4):492–500 Schweiger M, Arridge SR, Hiraoka M, Delpy DT (1995) The finite-element method for the propagation of light in scattering media—boundary and source conditions. Med Phys 22(11):1779– 1792 Srinivasan S, Pogue B, Dehghani H, Jiang S, Song X, Paulsen K (2004) Improved quantification of small objects in near-infrared diffuse optical tomography. J Biomed Opt 9(6):1161–1171 Tikhonov AN, Arsenin VY (1977) Solutions of ill-posed problems. Wiley, Washington, DC
