Improved sparse reconstruction for fluorescence molecular tomography with L 1 / 2 regularization

Fluorescence molecular tomography (FMT) is a promising imaging technique that allows in vivo visualization of molecular-level events associated with disease progression and treatment response. Accurate and efficient 3D reconstruction algorithms will facilitate the wide-use of FMT in preclinical research. Here, we utilize L1/2-norm regularization for improving FMT reconstruction. To efficiently solve the nonconvex L1/2norm penalized problem, we transform it into a weighted L1-norm minimization problem and employ a homotopy-based iterative reweighting algorithm to recover small fluorescent targets. Both simulations on heterogeneous mouse model and in vivo experiments demonstrated that the proposed L1/2-norm method outperformed the comparative L1-norm reconstruction methods in terms of location accuracy, spatial resolution and quantitation of fluorescent yield. Furthermore, simulation analysis showed the robustness of the proposed method, under different levels of measurement noise and number of excitation sources. ©2015 Optical Society of America OCIS codes: (170.6960) Tomography; (170.3010) Image reconstruction techniques; (170.3880) Medical and biological imaging; (100.3190) Inverse problems. References and links 1. A. Ale, V. Ermolayev, E. Herzog, C. Cohrs, M. H. de Angelis, and V. Ntziachristos, “FMT-XCT: in vivo animal studies with hybrid fluorescence molecular tomography-X-ray computed tomography,” Nat. Methods 9(6), 615– 620 (2012). 2. 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). 3. X. Montet, V. Ntziachristos, J. Grimm, and R. Weissleder, “Tomographic fluorescence mapping of tumor targets,” Cancer Res. 65(14), 6330–6336 (2005). 4. V. Ntziachristos, E. A. Schellenberger, J. Ripoll, D. Yessayan, E. Graves, A. Bogdanov, Jr., L. Josephson, and R. Weissleder, “Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin V-Cy5.5 conjugate,” Proc. Natl. Acad. Sci. U.S.A. 101(33), 12294–12299 (2004). 5. K. O. Vasquez, C. Casavant, and J. D. Peterson, “Quantitative whole body biodistribution of fluorescent-labeled agents by non-invasive tomographic imaging,” PLoS ONE 6(6), e20594 (2011). 6. S. Kossodo, M. Pickarski, S. A. Lin, A. Gleason, R. Gaspar, C. Buono, G. Ho, A. Blusztajn, G. Cuneo, J. Zhang, J. Jensen, R. Hargreaves, P. Coleman, G. Hartman, M. Rajopadhye, T. Duong, C. Sur, W. Yared, J. Peterson, and B. Bednar, “Dual in vivo quantification of integrin-targeted and protease-activated agents in cancer using fluorescence molecular tomography (FMT),” Mol. Imaging Biol. 12(5), 488–499 (2010). 7. K. Radrich, P. Mohajerani, J. Bussemer, M. Schwaiger, A. J. Beer, and V. Ntziachristos, “Limited-projectionangle hybrid fluorescence molecular tomography of multiple molecules,” J. Biomed. Opt. 19(4), 046016 (2014). 8. Y. Lin, W. C. Barber, J. S. Iwanczyk, W. W. Roeck, O. Nalcioglu, and G. Gulsen, “Quantitative fluorescence tomography using a trimodality system: in vivo validation,” J. Biomed. Opt. 15(4), 040503 (2010). 9. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007). 10. L. Zhao, H. Yang, W. Cong, G. Wang, and X. Intes, “Lp regularization for early gate fluorescence molecular tomography,” Opt. Lett. 39(14), 4156–4159 (2014). 11. G. Zhan, X. Cao, B. Zhang, F. Liu, J. Luo, and J. Bai, “MAP estimation with structural priors for fluorescence molecular tomography,” Phys. Med. Biol. 58(2), 351–372 (2013). 12. M. Solomon, R. E. Nothdruft, W. Akers, W. B. Edwards, K. Liang, B. Xu, G. P. Suddlow, H. Deghani, Y. C. Tai, A. T. Eggebrecht, S. Achilefu, and J. P. Culver, “Multimodal fluorescence-mediated tomography and SPECT/CT for small-animal imaging,” J. Nucl. Med. 54(4), 639–646 (2013). #234854 $15.00 USD Received 19 Feb 2015; revised 4 Apr 2015; accepted 5 Apr 2015; published 9 Apr 2015 (C) 2015 OSA 1 May 2015 | Vol. 6, No. 5 | DOI:10.1364/BOE.6.001648 | BIOMEDICAL OPTICS EXPRESS 1648 13. H. W. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems (Kluwer Academic Publishers, 2000). 14. 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), 3025–3035 (2009). 15. M. Freiberger, C. Clason, and H. Scharfetter, “Total variation regularization for nonlinear fluorescence tomography with an augmented Lagrangian splitting approach,” Appl. Opt. 49(19), 3741–3747 (2010). 16. J. Dutta, S. Ahn, C. Li, S. R. Cherry, and R. M. Leahy, “Joint L1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol. 57(6), 1459–1476 (2012). 17. A. Behrooz, H. M. Zhou, A. A. Eftekhar, and A. Adibi, “Total variation regularization for 3D reconstruction in fluorescence tomography: experimental phantom studies,” Appl. Opt. 51(34), 8216–8227 (2012). 18. J. Chamorro-Servent, J. F. Abascal, J. Aguirre, S. Arridge, T. Correia, J. Ripoll, M. Desco, and J. J. Vaquero, “Use of Split Bregman denoising for iterative reconstruction in fluorescence diffuse optical tomography,” J. Biomed. Opt. 18(7), 076016 (2013). 19. H. Yi, D. Chen, X. Qu, K. Peng, X. Chen, Y. Zhou, J. Tian, and J. Liang, “Multilevel, hybrid regularization method for reconstruction of fluorescent molecular tomography,” Appl. Opt. 51(7), 975–986 (2012). 20. Z. Xue, X. Ma, Q. Zhang, P. Wu, X. Yang, and J. Tian, “Adaptive regularized method based on homotopy for sparse fluorescence tomography,” Appl. Opt. 52(11), 2374–2384 (2013). 21. D. Han, J. Tian, S. Zhu, J. Feng, C. Qin, B. Zhang, and X. Yang, “A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization,” Opt. Express 18(8), 8630–8646 (2010). 22. J. C. Baritaux, K. Hassler, and M. Unser, “An efficient numerical method for general Lp regularization in fluorescence molecular tomography,” IEEE Trans. Med. Imaging 29(4), 1075–1087 (2010). 23. H. Yi, D. Chen, W. Li, S. Zhu, X. Wang, J. Liang, and J. Tian, “Reconstruction algorithms based on l1-norm and l2-norm for two imaging models of fluorescence molecular tomography: a comparative study,” J. Biomed. Opt. 18(5), 056013 (2013). 24. J. Zeng, J. Fang, and Z. Xu, “Sparse SAR imaging based on L1/2 regularization,” Sci. China. Inf. Sci. 55(8), 1755–1775 (2012). 25. R. Chartrand, “Exact reconstruction of sparse signals via nonconvex minimization,” IEEE Signal Process. Lett. 14(10), 707–710 (2007). 26. G. Gasso, A. Rakotomamonjy, and S. Canu, “Recovering sparse signals with a certain family of nonconvex penalties and DC programming,” IEEE Trans. Signal Process. 57(12), 4686–4698 (2009). 27. S. Okawa, Y. Hoshi, and Y. Yamada, “Improvement of image quality of time-domain diffuse optical tomography with l sparsity regularization,” Biomed. Opt. Express 2(12), 3334–3348 (2011). 28. X. Chen, D. Yang, Q. Zhang, and J. Liang, “L1/2 regularization based numerical method for effective reconstruction of bioluminescence tomography,” J. Appl. Phys. 115(18), 184702 (2014). 29. D. Zhu and C. Li, “Nonconvex regularizations in fluorescence molecular tomography for sparsity enhancement,” Phys. Med. Biol. 59(12), 2901–2912 (2014). 30. D. Zhu, Y. Zhao, R. Baikejiang, Z. Yuan, and C. Li, “Comparison of Regularization Methods in Fluorescence Molecular Tomography,” Photonics 1(2), 95–109 (2014). 31. Z. Xu, H. Guo, Y. Wang, and H. Zhang, “Representative of L1/2 Regularization among Lq (0<q<1) Regularizations: an Experimental Study Based on Phase Diagram,” Acta Automat. Sinica 38(7), 1225–1228 (2012). 32. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). 33. A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005). 34. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). 35. D. Wang, X. Liu, Y. Chen, and J. Bai, “A novel finite-element-based algorithm for fluorescence molecular tomography of heterogeneous media,” IEEE Trans. Inf. Technol. Biomed. 13(5), 766–773 (2009). 36. M. Asif and J. Romberg, “Fast and accurate algorithms for re-weighted L1-norm minimization,” IEEE Trans. Signal Process. 61(23), 5905–5916 (2013). 37. H. Zhang, Y. Wang, X. Chang, and Z. Xu, “L1/2 regularization,” Sci. China Ser. E 40(3), 412–422 (2010). 38. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007). 39. H. Guo, Y. Hou, X. He, J. Yu, J. Cheng, and X. Pu, “Adaptive hp finite element method for fluorescence molecular tomography with simplified spherical harmonics approximation,” J. Innov. Opt. Heal. Sci. 7(2), 1350057 (2014). 40. J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Process. 16(12), 2992–3004 (2007). 41. T. G. Feeman, The Mathematics of Medical Imaging (Springer New York, 2010), Chap. 9. 42. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,” Opt. Express 14(16), 7109–7124 (2006). 43. M. A. Naser and M. S. Patterson, “Improved bioluminescence and fluorescence reconstruction algorithms using diffuse optical tomography, normalized data, and optimized selection of the permissible source region,” Biomed. Opt. Express 2(1), 169–184 (2011). #234854 $15.00


Introduction
As a promising non-invasive molecular imaging technique, in vivo fluorescence molecular tomography (FMT) is an active research topic and exhibits significant potential in many preclinical research, such as monitoring enzyme activity [1], gene expression visualization [2], mapping expressions of cancer markers [3,4], and monitoring targeted drug delivery [5].FMT allows for three-dimensional localization and quantification of fluorescence distribution of targets in small animals to resolve biological processes at molecular and cellular levels [6][7][8][9].This takes place through modeling of near infrared (NIR) light propagation in biological tissues with a forward transportation model, and solving an inverse problem for recovering the distribution of fluorophore concentration from boundary measurements.
Mathematically, the inverse problem of FMT is inherently ill-posed due to more unknowns than number of detectors which makes FMT reconstruction subject to the effects of inevitable noise and numerical errors [10].The strong light absorption and scattering in biological tissues also add to the computational complexity.To obtain meaningful results researchers tend to incorporate some priori information such as anatomical information, optical properties of tissues or permissible region in FMT inverse problem [11,12].Moreover, various regularization techniques are essential to produce approximate stable reconstruction result [13].Conventionally, L 2 -norm regularization based reconstruction methods are applied to FMT.Although these methods can deal with the ill-posedness of FMT inverse problem, the over-smoothness of L 2 -norm results in blurred or spread targets with the loss of high-frequency feature in the reconstructed images [14].Owing to the edge-preserving properties of total variation (TV) norm, TV regularization based reconstruction methods are presented for linear and nonlinear FMT, and the numerical and experimental results therein demonstrate potential of offering higher resolution and robustness compared to conventional L 2 regularization [15][16][17].
Note that fluorescent probes in typical FMT applications are locally concentrated within specific regions of interest (e.g., inside tumors), sparse solutions are more preferred.Although L 0 -norm is an ideal sparsity regularizer, reconstruction based on L 0 -norm regularization involves a problem of combinatory optimization, which makes it inefficient or unfeasible in practical applications.Inspired by the theories of compressive sensing and sparse signal recovery, researchers resort to L 1 -norm based reconstruction methods and produce better results over L 2 -norm based methods in scenarios of recovering small localized targets [18][19][20][21][22][23].However, solutions of L 1 -norm methods are still less sparse and inefficient when heavytail distribution errors exist in the measurements [24].Recently, non-convex L p (0<p<1) norm regularizations have attracted much attention in the field of statistical learning and compressive sensing due to the ability of enhancing sparsity [25,26].Compared to L q (1≤q≤2) norm, L p (0<p<1) pseudo norm, as p→0 will be a better approximation to L 0 norm.Consequently, reconstruction methods with non-convex regularization are introduced into time-domain diffuse optical tomography (DOT) [27], bioluminescence tomography (BLT) [28], and time-domain FMT [10], providing improved reconstructions over L 1 -norm methods.Nonconvex regularizations for FMT are also investigated with numerical simulations and phantom experiments, where a majorization-minimization algorithm is adopted in the iterative reconstruction process in [29,30].Comparison studies in [30] further show nonconvex L p methods with p near 1/2 performs the best in reconstructing small targets, which is consistent with the conclusion in [31], i.e., L 1/2 regularizer is a good representative of L p (0<p<1) regularizers.
In this paper, we adopted an L 1/2 -norm regularization framework for FMT to recover the distribution of small fluorescent targets.To solve the nonconvex optimization problem introduced by L 1/2 regularizer, we convert it into a series of weighted L 1 -norm optimization problem and iteratively solve it by a homotopy-based algorithm.The performance of the proposed method in FMT reconstruction was validated with simulated data and in vivo experimental data.

Photon propagation model
In the NIR spectrum of 700-900 nm, diffusion equation is usually used for depicting photon propagation model in highly scattering media [32,33].For steady-state FMT with point excitation sources, the following coupled diffusion equations have been extensively used for forward model [33][34][35]: where , x m Φ denotes the photon flux density, and subscripts x and m denote the excitation and emission light, respectively., x m ax am sx sm D g ( ) af r ημ denotes the unknown fluorescent yield distribution to be reconstructed.l r represents the location of the point sources with an amplitude of s Θ .To solve these equations, the Robin-type boundary conditions are adopted on the boundary ∂Ω of the domain Ω [34].By solving the diffusion equation with the finite element method, we can obtain a linear model: where X is an 1 N × vector representing the unknown fluorescent source distribution, Φ is an 1 M × vector which contains the boundary measurements, and A is an M N × weight matrix, that depends on the geometry and the optical properties of the turbid medium.

The framework of L 1/2 -norm regularization
As mentioned above, a proper regularization is inevitable to deal with the high ill-posedness of FMT.With L p -norm regularization, an approximate solution to Eq. ( 2) can be obtained by solving the following minimization problem: where λ is the regularization parameter and P P X is the L p -norm of the solution X .By substituting p with 1/2, we get an objective function with L 1/2 -norm penalty term,  (1 / ) , where i x denotes the absolute value of i x .So Eq. ( 4) was converted into the following form: where 0 i w > denotes the weight at index i , and t x denotes the estimated signal in last iteration.To guarantee the feasibility, we used i for weights calculation, here κ is a small positive real number (e.g.1e 5 κ = − used in our implementation).Thus, by iterative reweighting, i.e. re-compute weights at every iteration using the solution at the previous iteration, we can successfully convert the non-convex L 1/2 -norm minimization problem into a series of weighted L 1 minimization problem, which can be solved by efficient L 1 minimization algorithms.Here, we adopted a homotopy-based algorithm for iterative reweighting that quickly updates the solution of Eq. ( 6) as the weights change and iteratively recover the unknown fluorescent distribution.Suppose that the weights update from t w to 1 t w + at the ( 1) t + th iteration.To incorporate changes in the weight and quickly compute the new solution of X , we propose the following homotopy program: where ε denotes the homotopy parameter.Obviously, by changing ε from 0 to 1, we can phase in the new weights and phase out the old ones.According to optimal theory, for any value of ε , the solution * X must satisfy the following optimality conditions [36]: ) for all where Γ denotes the support of * X , z denotes the sign sequence of * X of Γ , and i a denotes i th column of A .According to optimal theory, as we change ε to ε δ + , for a small value of δ , the solution moves in a direction X ∂ and the optimality conditions change as [36]: where t w and 1 t w + denote Γ × Γ diagonal matrices with their diagonal entries being the values of t w and With the change of ε , by the new optimality conditions Eq. (9.a) the update direction of As we increase δ , the solution moves in the direction X ∂ until either a new element enters the support of the solution (when a constraint in Eq. (9.b) becomes active), or an existing element shrinks to zero.The step size that takes the solution to such a critical value of ε can be computed as * min( , ) δ δ δ + − = where: * min( , ) , min( ) δ + denotes the smallest step-size that causes a constraint in Eq. (9.b) to become active, indicating entry of a new element at index γ + in the support, whereas δ − denotes the smallest step-size that shrinks an existing element at index γ − to zero.
The new critical value of ε becomes ε δ , where its support and sign sequence are updated accordingly.At every homotopy step, we jump from one critical value of ε to the next while updating the support of the solution, until ε = 1.
In our implementation, we set the initial solution as 0 (1, ,1) X =  , which means the solution of the first iteration corresponds to a solution of L 1 -norm minimization.The iterative reconstruction algorithm stops when the number of iteration reaches a given value k.

Experiments and results
In this section, both simulated data on a 3D digital mouse model and in vivo experimental data acquired by a dual-modality FMT/Micro-CT system were used to validate the potential and feasibility of the proposed iterative reweighted homotopy based L 1/2 -norm regularization algorithm, which we will call IRW-L 1/2 .Several metrics were used to evaluate the quality of the reconstructed results, i.e. the location error (LE), quantification error (QE), and the percentage of non-zero coefficient (PNZ).LE is the Euclidean distance between the center of the reconstructed and the actual target.QE is defined as the relative error between the reconstructed fluorescent yield (RFY) and the true value.PNZ is the percentage of non-zero components in the vector of solution, which can be regarded as a metric of sparse degree of the solution.
A series of verification experiments were designed in section 3. First, a single-target reconstruction was performed to verify the reconstruction accuracy of our method in target location and fluorescent yield.Then, the influence of target depth and relative position of the excitation nodes on the performance of reconstruction algorithms was thoroughly investigated by a comparative study in section 3.2.In section 3.3, we analyzed the robustness and stability of the proposed algorithm over noise or reduction of excitation sources.The multiple-target resolving ability of the proposed method was also evaluated with double-target experiments in section 3.4.Finally, an in vivo experiment was conducted to further test our method.The reconstruction code written in MATLAB was performed on our desktop computer (Intel ® Xeon ® Processor E3-1230).In our implementation, the regularization parameters for all the tested L 1/2 -norm and L 1 -norm methods were selected empirically, ranging from 1e-9 to 1e-14.The maximum iteration number k = 50.

Single-target reconstruction on 3D digital mouse atlas
In the numerical simulations, a 3D digital mouse atlas was utilized to provide anatomical information [38].Only the torso with a height of 40mm was investigated in our experiments, which consisted of six organs, i.e., muscle, heart, liver, stomach, kidneys and lungs.Table 1 lists the related optical properties [39].
A cylindrical fluorescent target with a 0.8mm radius and 1.6mm height was placed in the liver with center at P 0 = (11.9,6.4, 16.4mm), as shown in Fig. 1(a).The fluorescent yield of the target was set to be 1 0.05mm − .Figure 1(c) shows the distribution of the 18 point sources used in this simulation.The black dots represent positions of the point excitation sources, which were located one transport mean free path beneath the surface on the plane of Z = 16.4mm.For each excitation source, the surface data on the opposite side with a 120° field of view (FOV) were measurable.Measurements were obtained every 20° and a total of 18 data sets were assembled for the reconstruction of the fluorescent yield.For the forward calculation, the torso was discretized into a FEM mesh with 19278 nodes and 101329 tetrahedral elements.Figure 1(b) shows the simulated measurements on the surface.For the inverse problem, a different mesh with 3051 nodes and 16453 tetrahedral elements was used for FMT reconstruction.We thoroughly compared six methods in single-target reconstruction, including the proposed IRW-L 1/2 , L 1 -norm based homotopy algorithm (Homotopy-L 1 ), the incomplete variables truncated conjugate gradient algorithm (IVTCG) [23], TV-norm based two-step iterative shrinkage algorithm [40], L 2 -norm based Tikhonov regularization, and the algebraic reconstruction technique (ART) [41].Figure 2 shows the cross-section views and 3D renderings of the reconstruction by six methods.Table 2 summarizes the detailed results.Compared to the other four methods, the reconstructed images by IRW-L 1/2 and Homotopy-L 1 have fewer artifacts by visual assessment, and the corresponding location errors are lower.Furthermore, we found that the PNZ corresponding to IRW-L 1/2 is the lowest, which means the solution by our proposed method is the sparsest.As for the reconstructed fluorescent yield, the L 1/2 -norm method outperforms other methods.In general, the proposed IRW-L 1/2 method produced the best overall result in terms of LE, QE and PNZ in this case.It is worth mentioning that we did not combine any auxiliary strategy, such as permissible regions or multilevel FEM, with the six methods to improve reconstruction.

Influence of depth and relative position of the excitation nodes on algorithm performance
We also investigated the variation of reconstruction accuracy with target depth and relative position of the excitation node.As shown in Fig. 3, the following cases were considered: (a) keep the excitation sources fixed and move the single-target along Y-axis direction from P0 to different positons P1, P2 P3, and P4; (b) Fix the target position at P3 and move the excitation sources from plane of Z0 to planes of Z1, Z2 and Z3.We employed the above six methods to reconstruct the single-target at four different positions in case (a) and displayed the comparison results with 3 histograms.As shown in Fig. 4(a) and Fig. 4(b), IRW-L 1/2 produced the lowest location error and highest fluorescent yield in all of the test situations.We see the location error by IRW-L 1/2 method is within 0.65mm at different positions.It indicates the location accuracy by IRW-L 1/2 method varies slightly with target depth.Here, we also calculated the value of full width at half maximum (FWHM) to measure the retrieved spatial resolution, as shown in Fig. 4(c).We found that results corresponding to IRW-L 1/2 and Homotopy-L 1 have lower FWHM, which means they produce improved spatial resolution over the other methods.
Figure 5 compared the results of the six algorithms for reconstructing the single-target at P3 with the excitation sources located at different planes.We found that the relative position of excitation nodes to target has more influence on the retrieved fluorescent yield and location accuracy than on FWHM.We see an obvious increase of LE from 0.1mm to 0.91mm when the excitation sources moved from plane of Z0 to planes of Z1, Z2 and Z3.However, the overall comparative results in Fig. 4 and Fig. 5 reveal that IRW-L 1/2 performs best and Homotopy-L 1 takes the second place in most of the test cases.Therefore, only IRW-L 1/2 and Homotopy-L 1 are considered in the following simulations.

Stability analysis
To evaluate the robustness and stability of the proposed reconstruction algorithm, we took into consideration the influence of noise and the number of excitation sources on FMT reconstruction in this section.
Generally speaking, the quality of FMT reconstruction is sensitive to noise in measurements due to the ill-posedness of inverse problem.By adding different levels (0%-30%) of Gaussian noise to the simulated data, we investigated the performance of the proposed IRW-L 1/2 reconstruction algorithm on noisy data.Figure 6 illustrates variations of LE and the absolute error in the reconstructed fluorescent yield at different noise levels.Although both the results of Homotopy-L 1 and IRW-L 1/2 algorithm suffer from some degree of fluctuation, the performance of IRW-L 1/2 is more stable than L 1 -norm with the increase of noise level.Especially, the LE by IRW-L 1/2 in all of the testing cases is within 0.7mm, which is smaller than that of Homotopy-L 1 .On the other side, reduction in measurements that result from decreasing excitation sources will also aggravate the ill-posed inverse problem and thus affect the reconstruction.In section 3.1 and 3.2, excitation nodes were obtained every 20° and a total of 18 data sets were assembled for the subsequent reconstruction process.In this section, by gradually and uniformly reducing excitation nodes from 18 to 1, we investigated the dependence of the reconstruction accuracy on the number of excitation sources.The illustration in Fig. 7 shows the comparison results.Obviously, the fluctuation caused by the reduction of excitation nodes looks bigger than by noise in measurements, especially for L 1 -norm method.However, IRW-L 1/2 is less sensitive to reduction of measurements and outperforms Homotopy-L 1 .Even in the case of only one excited node, our method still produces a satisfactory result with a LE of 0.81mm.In contrast, when the excitation node number goes smaller than 6, performance of Homotopy-L 1 deteriorates rapidly, with location deviation more than 1mm.
The stability analysis in this section demonstrates that the performance of the proposed reconstruction algorithm IRW-L 1/2 is superior to Homotopy-L 1 .As shown in Fig. 6 and Fig. 7, the proposed method produced very stable results from noisy or limited measurements with high location accuracy within 0.81mm.

Reconstruction of double-target
Accurate resolving multiple targets is difficult and crucial for FMT reconstruction algorithms.In this section, we evaluate the multiple-target resolving ability of the proposed algorithm with two groups of double-target experiments on the same mouse model in section 3.1.
In the first group test, two identical targets were placed at T1 = (11.8,6.4, 16.4mm) and T2 = (11.8,10.9, 16.4mm), respectively, with a center-to-center separation of 4.5mm, as shown in Fig. 8(d).The size of the targets is same as that in single-target case.Figure 8(a) illustrates the simulated photon distribution on the forward mesh.For reconstruction, the mouse model was discretized into a FEM mesh with 3960 nodes and 21539 tetrahedral elements.Figure 8 presents the related results by the two reconstruction algorithms.It is clear that the solution of IRW-L 1/2 is concentrated in two localized regions and with two targets clearly separated, as depicted in Fig. 8(b) and Fig. 8(e).The detailed comparison results in Table 3 demonstrate that IRW-L 1/2 outperform the Homotopy-L 1 in terms of all of the metrics.Meanwhile, the reconstructed two targets are very similar in size and fluorescence yield, and the location errors for two targets are 0.78 and 1.02mm, respectively.In contrast, although two targets could be resolved by Homotopy-L 1 , the related result has larger deviation in target positions and fluorescence yield.In the second group test, we considered a more challenging situation, i.e. two closer targets were placed in the liver at T1 = (11.8,10.8, 16.4mm) and T2 = (11.8,13.8, 16.4mm), respectively.The distance between two targets is 3mm, as shown in Fig. 9(d).Considering that the radius of the targets is 0.8mm, the actual boundary-to-boundary distance between T1 and T2 is only 1.4mm.In this case, T2 is a position close to back, the simulated photon distribution on the surface is quite different from that in the above case (separation of 4.5mm), as shown in Fig. 8(a) and Fig. 9 4 summarizes the detailed quantitative results in the second group test.As we can see in Fig. 9, our method successfully recover two targets, while Homotopy-L 1 fail to discriminate such close two targets.It means that IRW-L 1/2 provides better spatial discrimination than Homotopy-L 1 .It is worth explaining the reconstructed center in Table 3 and Table 4.For better evaluation the location accuracy of reconstruction algorithm, we need to resolve the individual center of targets.In single-target case, we took the nodes with maximum value as the center position of the target.Nevertheless, we cannot deal with the solution in this direct way in double-target cases.To resolve the solution, the nodes with value bigger than 10% of the maximum reconstructed value were chosen, and then two clusters were constructed according to the coordinates of these nodes (using a MATLAB function clusterdata).For each cluster, we took the nodes with the maximum reconstructed value as one reconstructed center.

In vivo evaluation with implanted fluorophore
In this section, we further assess the performance of the proposed algorithm with in vivo experimental data, which come from [23].The in vivo experiment was conducted on an adult BALB/C mouse with a dual-modality FMT/Micro-CT system.In this experiment, a glass tube with 0.6 mm radius and 2.8 mm height, which contains Cy5.5 solution (with the extinction coefficient of about 0.019mm -1 µM -1 and quantum efficiency of 0.23 at the peak excitation wavelength of 671 nm) [42], was implanted into the abdomen of the anesthetized mouse.The true fluorescent yield of Cy5.5 is 1 0.0402mm − [43].With Micro-CT, the true center of the target was determined as (28.5, 25.7, 13.4mm).The CT data was segmented into five major anatomical components, including heart, lungs, liver, kidneys, and muscle, as shown in Fig. 10(e).Table 5 lists the optical properties at both the excitation and emission wavelengths, which were calculated based on the literature [44].Figure 10(a) shows the measurements mapped onto the 3D surface of the mouse torso.For subsequent reconstruction, the segmented mouse torso was discretized into a mesh with 18,504 tetrahedral elements and 3823 nodes. and

Fig. 1 .
Fig. 1.(a) Torso of the mouse atlas model with a cylindrical fluorescent targets in the liver, (b) Forward mesh and the simulated photon distribution on surface.(c) Point excitation sources on the plane of 16.4 Z mm = .The 18 black points denotes the point excitation sources.For each excitation source, fluorescence data is detected at the opposite side with a 120° FOV.

Fig. 2 .
Fig. 2. Reconstruction results in single-target case.(a) Cross-section views of the reconstructions by six comparative methods at the axial slice where the center of the real fluorescent target (denoted by the small black circle) is located; (b) 3D renderings of the results.The red cylinder is the real target while the green zone denotes the reconstructed target.

Fig. 4 .
Fig. 4. Comparison results for reconstruction single-target at different depths, in terms of (a) location error, (b) fluorescent yield, and (c) full width at half maximum.

Fig. 5 .
Fig. 5. Comparison result for reconstruction a single-target with the excitation sources located at different planes, in terms of (a) location error, (b) fluorescent yield, and (c) full width at half maximum.

Fig. 6 .
Fig. 6.Illustrations of variations of error of RFY (a) and LE (b) obtained at different noise levels.Subfigure (a) and (b) correspond to absolute quantitative error of reconstructed fluorescent yield and location error, respectively.

Fig. 7 .
Fig. 7. Illustration of variations of R.FY (a) and LE (b) with different number of excitation nodes.Subfigures (a) and (b) correspond to absolute quantitative error of reconstructed fluorescent yield and location error, respectively.
(a).For inverse calculation, the digital mouse model was discretized into a mesh with 5951 nodes and 29671 tetrahedral elements.The final reconstruction by the IRW-L 1/2 is shown in Fig. 9(b) and Fig. 9(e).For comparison, Fig. 9(c) and Fig. 9(f) present the result by Homotopy-L 1 .Table

Fig. 8 .
Fig. 8. Reconstruction results of two comparative methods for double targets with 4.5mm separation.Subfigure (a) shows the forward mesh and the photon distribution; (d) shows the digital mouse model with two targets with 4.5mm separation; (b) and (e) correspond to the results of IRW-L 1/2 , (c) and (f) present the results of Homotopy-L 1 ; (b) and (c) correspond to the transverse view at the plane of z = 16.4mm,where the small black circle denotes the true source; (e) and (f) are the corresponding 3D view, where the red cylinder denotes the true target while the green zone is the reconstructed target.

Fig. 9 .
Fig. 9. Reconstruction results of two comparative methods for double fluorescent targets with 3mm separation.Subfigure (a) shows the forward mesh and the photon distribution; (d) shows the digital mouse model with two targets with 3mm separation; (b) and (e) correspond to the results of IRW-L 1/2 ; (c) and (f) present the results of Homotopy-L 1 .(b) and (c) correspond to the transverse view at the plane of z = 16.4mm,where the small black circle denotes the actual source; (e) and (f) illustrate the 3D view of the reconstruction, where the red cylinder denotes the true target while the green zone is the reconstructed target.