1.
Introduction
From being the solid development of bioluminescent probes and reporter technologies,^{1, 2, 3} the application of bioluminescence in biomedical in vivo imaging has become more and more attractive over the recent years. It offers an alternative opportunity for noninvasively visualizing biological processes at the physiological and molecular levels in whole animals.^{4, 5, 6} Tomographic bioluminescence imaging (TBI) can further translate the planar imaging information into three-dimensional bioluminescent source distribution quantitatively, thus greatly facilitating applications in related biomedical in vivo studies.^{7, 8} Therefore, the development of the reconstruction approaches for tomographic imaging plays a unique role in the achievement of practical biomedical in vivo imaging studies.
Nevertheless, it is known that the inverse problem of tomographic imaging is an ill-posed problem due to the fact that some limited information can only be measured from the boundary of animals to estimate the internal bioluminescent source distribution.^{9} Therefore, this imaging modality faces various challenges in its accuracy and reliability. Researchers have made considerable efforts to cope with the ill-posedness of the underdetermined inverse problem. It can be alleviated by incorporating some a priori information. Spectrally resolved boundary measurements, which practically increase the amount of independent data, are necessary to accurately recover tomographic images of bioluminescent sources.^{10, 11, 12, 13, 14} By reducing the number of unknowns, the permissible source region method is demonstrated to be capable of correctly reconstructing the images at monochromatic measurements.^{15, 16, 17} This reconstruction domain is needed to be spatially constrained to the area of interest, whereas, it is not always feasible to define such a region effectively in practical tomographic imaging applications. Moreover, the varying boundary conditions have also been proposed to enhance the reliability of the reconstruction results.^{18}
On the other hand, the regularization strategies, which are imposed with the output-least-squares formulation to stabilize the inverse problem, are also indispensable for tomographic quality.^{19} Although the l _{2}-type regularization strategy is the most popular and commonly-applied, it often imposes over-smoothing on the reconstruction results. In contrast, due to the high specificity of the bioluminescence probes and reporter technologies, the objective for tomographic imaging is only the sparsely-distributed signal. Consequently, the sparse prior knowledge can be employed to recover the source distribution and preserve discontinuities in the reconstructed profiles with less measurements, which is presented in references by Lu ^{20} and Gao and Zhao.^{21} Actually, the sparse-type (l _{1}) regularization has been studied for years and recently has drawn a lot of attention due to theoretical justification and many applications. ^{20, 22, 23, 24, 25, 26, 27} Although the aforementioned methods have been proposed to overcome the challenges in the inverse problem, further study is still needed to alleviate the dependence with some parameters during reconstructions, such as the parameters in the regularization term.
In general, it is assumed that there are multiple local minima of the objective function in the TBI inverse problem.^{28} Consequently, the application of local optimization techniques, e.g., the conventional Newton-based optimization methods that explore the unknown parameter space only near a single local minimum, may lead to inapposite convergence without finding the global minimum.^{16, 17, 29, 30} As a consequence, reconstructions on a part of the animal body or even on the whole body become more difficult. Recently, a graph cuts-based method was proposed, which could find a global optimal solution efficiently and was not dependent on a starting unknown guess in the search process.^{31} However, since the reconstruction could only be represented in discrete form, it was difficult to recover complex source distribution appropriately in whole animals. Likewise, more efforts are urgently needed to make whole body imaging available with an arbitrary initial guess for the unknowns.
In this work, a dynamically-sparse regularized global method is proposed. In order to make full use of the sparse a priori information, the sparse term l _{p} (1 ⩽ p < 2) is approximated by a corresponding weighted l _{2} norm in each iteration,^{32} rather than be adopted directly. It can facilitate the generalization of the sparse regularization with greater flexibility of p in the dynamic quadratic frame, and avoid tedious numerical operation as well. A globally-converged optimization technique is presented to search for the global solution. It could find the globally optimal solutions far from the starting unknown guess efficiently, and maintain the reliable reconstruction quality over a wide-range of regularization parameters.
In Sec. 2, we present the proposed method for TBI. In Sec. 3, validations based on an in vivo experimental data set demonstrate higher imaging reliability and cheaper computational cost of the proposed method. First reconstruction comparisons of the proposed method with other methods demonstrate the applicability on the entire region. Then, the reliable performance on various ranges of regularization parameters and initial unknown values is also validated. Based on the in vivo mouse experiment and a mouse atlas, the tolerance for optical property mismatch is investigated with optical overestimation and underestimation on absorbing and reduced scattering properties. Additionally, the reconstruction efficiency is also analyzed with different sizes of mouse grids. Finally, we discuss and conclude this paper.
2.
Methods
2.1.
Forward Model
In the steady-state domain, the forward problem of light propagation for TBI can be modeled as a diffusion approximation to the radiative transport equation,^{16} which is given by
Eq. 1
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} -\nabla \cdot (D(\textbf {r})\nabla [\Phi (\textbf {r})]+ \mu _{a}(\textbf {r})\Phi (\textbf {r})=\mathcal {X}(\textbf {r}) (\textbf {r}\in \Omega). \end{equation} \end{document} $$-\nabla \xb7(D\left(\mathbf{r}\right)\nabla \left[\Phi \left(\mathbf{r}\right)\right]+{\mu}_{a}\left(\mathbf{r}\right)\Phi \left(\mathbf{r}\right)=\mathcal{X}\left(\mathbf{r}\right)(\mathbf{r}\in \Omega ).$$In the bioluminescence imaging experiments, the whole process is performed in a completely dark environment, and there is no external photon into imaging domain Ω through its boundary ∂Ω, so the Robin-type boundary condition is suitable:
Eq. 2
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \Phi (\textbf {r})+2\kappa (\textbf {r},n,n^{\prime })D(\textbf {r}) [\textbf {v}(\textbf {r})\cdot \nabla \Phi (\textbf {r})] =0 (\textbf {r}\in \partial \Omega). \end{equation} \end{document} $$\Phi \left(\mathbf{r}\right)+2\kappa (\mathbf{r},n,{n}^{\prime})D\left(\mathbf{r}\right)[\mathbf{v}\left(\mathbf{r}\right)\xb7\nabla \Phi \left(\mathbf{r}\right)]=0(\mathbf{r}\in \partial \Omega ).$$Eq. 3
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} V(\textbf {r})=-D(\textbf {r}) [\textbf {v}(\textbf {r})\cdot \nabla \Phi (\textbf {r})] =\frac{\Phi (\textbf {r})}{2\kappa (\textbf {r},n,n^{\prime })} (\textbf {r}\in \partial \Omega). \end{equation} \end{document} $$V\left(\mathbf{r}\right)=-D\left(\mathbf{r}\right)[\mathbf{v}\left(\mathbf{r}\right)\xb7\nabla \Phi \left(\mathbf{r}\right)]=\frac{\Phi \left(\mathbf{r}\right)}{2\kappa (\mathbf{r},n,{n}^{\prime})}(\mathbf{r}\in \partial \Omega ).$$Based on the finite element theory, Eqs. 1, 2 are discretized by piecewise linear finite elements, and a matrix-vector equation is integrally assembled.^{33} Consequently, after a series of rearrangements for the elements in the matrix, the linear relationship between the measured photon density distribution and the unknown source distribution in heterogeneous biological tissues is established as:
Eq. 4
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \mathcal {MX}=\Phi, \end{equation} \end{document} $$\mathcal{MX}=\Phi ,$$2.2.
Dynamic Sparse Regularized Function
In practice, the common approach for the TBI inverse problem is to adopt the output-least-squares formulation. A regularization function is incorporated to stabilize the inversion problem like this:
Eq. 5
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} S(\mathcal {X})=||\mathcal {LX}||_{p}^{p}, \end{equation} \end{document} $$S\left(\mathcal{X}\right)={\left|\right|\mathcal{LX}\left|\right|}_{p}^{p},$$Eq. 6
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} T(\mathcal {X})=\frac{1}{2}\Vert \mathcal {MX}-\Phi ^{m}\Vert _{2}^{2}+\lambda ||\mathcal {LX}||_{p}^{p}, \end{equation} \end{document} $$T\left(\mathcal{X}\right)=\frac{1}{2}\parallel \mathcal{MX}-{\Phi}^{m}{\parallel}_{2}^{2}+{\lambda \left|\right|\mathcal{LX}\left|\right|}_{p}^{p},$$Eq. 7
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} S(\mathcal {X})=\frac{1}{p}\Vert \mathcal {X}\Vert _{p}^{p} (1\le p< 2). \end{equation} \end{document} $$S\left(\mathcal{X}\right)=\frac{1}{p}{\parallel \mathcal{X}\parallel}_{p}^{p}(1\le p<2).$$Now let us begin to construct the dynamic sparse regularized function. Instead of solving the sparse l _{p} problem,^{20} a quadratic version is generated to approximate the sparse term for each iteration.^{32} In order to replace the l _{p} term by the l _{2} one and dynamically regularize the objective function, the quadratic function is defined as:
Eq. 8
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} Q_{S}^{(k)}(\mathcal {X})=\frac{1}{2}\big \Vert W_{S}^{(k)^{1/2}}\mathcal {X}\big \Vert _{2}^{2} +\left (1-\frac{p}{2}\right)S(\mathcal {X}^{(k)}), \end{equation} \end{document} $${Q}_{S}^{\left(k\right)}\left(\mathcal{X}\right)=\frac{1}{2}\parallel {W}_{S}^{{\left(k\right)}^{1/2}}\mathcal{X}{\parallel}_{2}^{2}+\left(1-\frac{p}{2}\right)S\left({\mathcal{X}}^{\left(k\right)}\right),$$Eq. 9
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} W_{S}^{(k)}=\hbox{diag} [\tau _{S,\epsilon _{S}}(\mathcal {X}^{(k)})]. \end{equation} \end{document} $${W}_{S}^{\left(k\right)}=\text{diag}\left[{\tau}_{S,{\epsilon}_{S}}\left({\mathcal{X}}^{\left(k\right)}\right)\right].$$The diagonal matrix [TeX:] $W_{S}^{(k)}$ ${W}_{S}^{\left(k\right)}$ is updated in each iteration using the values of the last step. The matrix actually plays a spatially varying role to maintain the imaging results reliable enough, which can tolerate different orders of magnitude of regularization parameters to obtain reliable results to some extent. Based on the strategy in the reference by Rao and Kreutz-Delgato,^{37} [TeX:] $\tau _{S,\epsilon _{S}}$ ${\tau}_{S,{\epsilon}_{S}}$ (for some small ε_{S}) is defined as
Eq. 10
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \tau _{S,\epsilon _{S}}(u) = \left\lbrace \begin{array}{@{}l@{\quad}l} \mid u\mid ^{p-2} & \hbox{if\quad $\mid u \mid >\epsilon _{S}$} \\[3pt] 0 & \hbox{if\quad $\mid u \mid \le \epsilon _{S}$} \end{array}. \right. \end{equation} \end{document} $${\tau}_{S,{\epsilon}_{S}}\left(u\right)=\left\{\begin{array}{cc}\hfill \mid u{\mid}^{p-2}& \hfill \text{if}\phantom{\rule{1em}{0ex}}\mid u\mid >{\epsilon}_{S}\\ \hfill 0& \hfill \text{if}\phantom{\rule{1em}{0ex}}\mid u\mid \le {\epsilon}_{S}\end{array}.\right.$$Note that it is necessary to add the constant term [TeX:] $(1-\frac{p}{2})S(\mathcal {X}^{(k)})$ $(1-\frac{p}{2})S\left({\mathcal{X}}^{\left(k\right)}\right)$ in Eq. 8 to ensure that:
Eq. 11
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} S(\mathcal {X}^{(k)})=Q_{S}^{(k)}(\mathcal {X}^{(k)}) (\hbox{as} \epsilon _{S}\rightarrow 0), \end{equation} \end{document} $$S\left({\mathcal{X}}^{\left(k\right)}\right)={Q}_{S}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right)(\text{as}{\epsilon}_{S}\to 0),$$Eq. 12
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} S(\mathcal {X})<Q_{S}^{(k)}(\mathcal {X}) (\forall \mathcal {X}\ne \mathcal {X}^{(k)}), \end{equation} \end{document} $$S\left(\mathcal{X}\right)<{Q}_{S}^{\left(k\right)}\left(\mathcal{X}\right)(\forall \mathcal{X}\ne {\mathcal{X}}^{\left(k\right)}),$$Eq. 13
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \nabla _{\mathcal {X}} S(\mathcal {X})\mid _{\mathcal {X}=\mathcal {X}^{(k)}}=\nabla _{\mathcal {X}} Q_{S}^{(k)}(\mathcal {X})\mid _{\mathcal {X}=\mathcal {X}^{(k)}}\!\!. \end{equation} \end{document} $${\nabla}_{\mathcal{X}}S\left(\mathcal{X}\right){\mid}_{\mathcal{X}={\mathcal{X}}^{\left(k\right)}}={\nabla}_{\mathcal{X}}{Q}_{S}^{\left(k\right)}\left(\mathcal{X}\right){\mid}_{\mathcal{X}={\mathcal{X}}^{\left(k\right)}}.$$It is noteworthy that the original regularization term [Eq. 7] and its quadratic version [Eq. 8] have the same value and tangent direction at [TeX:] $\mathcal {X}=\mathcal {X}^{(k)}$ $\mathcal{X}={\mathcal{X}}^{\left(k\right)}$ .
Thus, based on the dynamic sparse regularization technique, the objective function is reformulated into the following quadratic form:
Eq. 14
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} T^{(k)}(\mathcal {X})&=&\frac{1}{2}\Vert \mathcal {MX}-\Phi ^{m}\Vert _{2}^{2}+\lambda Q_{S}^{(k)}(\mathcal {X})\nonumber \\ &=&\frac{1}{2}\Vert \mathcal {MX}-\Phi ^{m}\Vert _{2}^{2}+\frac{\lambda }{2}\big \Vert W_{S}^{(k)^{1/2}}\mathcal {X}\big \Vert _{2}^{2}\nonumber \\ &+&\lambda \left(1-\frac{p}{2}\right)S(\mathcal {X}^{(k)}) (k\ge 0). \end{eqnarray} \end{document} $$\begin{array}{ccc}\hfill {T}^{\left(k\right)}\left(\mathcal{X}\right)& =& \frac{1}{2}{\parallel \mathcal{MX}-{\Phi}^{m}\parallel}_{2}^{2}+\lambda {Q}_{S}^{\left(k\right)}\left(\mathcal{X}\right)\hfill \\ & =& \frac{1}{2}{\parallel \mathcal{MX}-{\Phi}^{m}\parallel}_{2}^{2}+\frac{\lambda}{2}\parallel {W}_{S}^{{\left(k\right)}^{1/2}}\mathcal{X}{\parallel}_{2}^{2}\hfill \\ & +& \lambda \left(1-\frac{p}{2}\right)S\left({\mathcal{X}}^{\left(k\right)}\right)(k\ge 0).\hfill \end{array}$$The corresponding gradient and Hessian matrix are easily derived
15.
Eq. 15a
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \nabla T^{(k)}(\mathcal {X})=\big(\mathcal {M}^{T}\mathcal {M}+\lambda W_{S}^{(k)}\big)\mathcal {X}-\mathcal {M}^{T}\Phi ^{m}, \end{equation} \end{document} $$\nabla {T}^{\left(k\right)}\left(\mathcal{X}\right)=\left({\mathcal{M}}^{T}\mathcal{M}+\lambda {W}_{S}^{\left(k\right)}\right)\mathcal{X}-{\mathcal{M}}^{T}{\Phi}^{m},$$Eq. 15b
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \nabla ^{2} T^{(k)}(\mathcal {X})=\mathcal {M}^{T}\mathcal {M}+\lambda W_{S}^{(k)}. \end{equation} \end{document} $${\nabla}^{2}{T}^{\left(k\right)}\left(\mathcal{X}\right)={\mathcal{M}}^{T}\mathcal{M}+\lambda {W}_{S}^{\left(k\right)}.$$In particular, if p = 2 and ε_{S} = 0, [TeX:] $Q_{S}^{(k)}(\mathcal {X})=\frac{1}{2}\Vert \mathcal {X} \Vert _{2}^{2}$ ${Q}_{S}^{\left(k\right)}\left(\mathcal{X}\right)=\frac{1}{2}{\parallel \mathcal{X}\parallel}_{2}^{2}$ for any k, thus the dynamically sparse regularized function degrades into the popular static regularized quadratic objective function.^{16, 20, 31}
2.3.
Reconstruction Algorithm
After constructing the dynamically sparse regularized function for the inverse problem, the reconstructed solution can be estimated by minimizing the objective function dynamically:
Eq. 16
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \mathcal {X}_{recons.}=\arg \mathop {\min _{\mathcal {X}}}T^{(k)}(\mathcal {X}), \end{equation} \end{document} $${\mathcal{X}}_{recons.}=\mathrm{arg}\underset{\mathcal{X}}{\mathrm{min}}{T}^{\left(k\right)}\left(\mathcal{X}\right),$$Eq. 17
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \nabla T^{(k)}(\mathcal {X})=0. \end{equation} \end{document} $$\nabla {T}^{\left(k\right)}\left(\mathcal{X}\right)=0.$$One kind of widely applied algorithm for solving Eq. 17 is the Newton method, as shown in Algorithm: where [TeX:] $\textbf {r}_{k}$ ${\mathbf{r}}_{k}$ denotes the increment step at kth iteration. This method exactly solves Newton's equations (step 3 in Algorithm) at each iteration, which can be very expensive if the number of unknowns is large and may not be justified when [TeX:] $\textbf {r}_{k}$ ${\mathbf{r}}_{k}$ is far from a solution. Thus, an inexact Newton method is preferred to just compute an approximate solution of Newton's equations at each iteration, which can be summarized in pseudo-code as follows:^{38}
Table 1
Newton method.
1: Initialize [TeX:] $\mathcal {X}^{(0)}$ ${\mathcal{X}}^{\left(0\right)}$ and tol |
2: while have not converged do |
3: Solve [TeX:] $\nabla ^{2} T^{(k)}(\mathcal {X}^{(k)})\textbf {r}_{k}=-\nabla T^{(k)}(\mathcal {X}^{(k)})$ ${\nabla}^{2}{T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right){\mathbf{r}}_{k}=-\nabla {T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right)$ |
4: Set [TeX:] $\mathcal {X}^{(k+1)}=\mathcal {X}^{(k)}+\textbf {r}_{k}$ ${\mathcal{X}}^{(k+1)}={\mathcal{X}}^{\left(k\right)}+{\mathbf{r}}_{k}$ |
5: end while |
Such approximate treatmental offers a trade-off between the accuracy with a solution of Newton's equations and the amount of the computational cost per iteration. By making the η_{k} appropriately small, the convergence can be made and the norm of [TeX:] $\nabla ^{2} T^{(k)}(\mathcal {X})$ ${\nabla}^{2}{T}^{\left(k\right)}\left(\mathcal{X}\right)$ can be reduced.^{39} Since [TeX:] $\nabla ^{2} T^{(k)}(\mathcal {X}^{(k)})\textbf {r}_{k}\break +\nabla T^{(k)}(\mathcal {X}^{(k)})$ ${\nabla}^{2}{T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right){\mathbf{r}}_{k}+\nabla {T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right)$ is the residual of Newton's equations, each η_{k} denotes how accurate [TeX:] $\textbf {r}_{k}$ ${\mathbf{r}}_{k}$ is close to the exact solution. It is seen that [TeX:] $\textbf {r}_{k}$ ${\mathbf{r}}_{k}$ satisfying step 4 in Algorithm is an inexact Newton step.
Table 2
Inexact Newton method.
1: Initialize [TeX:] $\mathcal {X}^{(0)}$ ${\mathcal{X}}^{\left(0\right)}$ and η_{0} |
2: while have not converged do |
3: Find some η_{k} ∈ [0, 1) and [TeX:] $\textbf {r}_{k}$ ${\mathbf{r}}_{k}$ that satisfy |
4: [TeX:] $\Vert \nabla ^{2} T^{(k)}(\mathcal {X}^{(k)})\textbf {r}_{k}+\nabla T^{(k)}(\mathcal {X}^{(k)})\Vert \le \eta _{k}\Vert \nabla T^{(k)}(\mathcal {X}^{(k)})\Vert$ $\parallel {\nabla}^{2}{T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right){\mathbf{r}}_{k}+\nabla {T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right)\parallel \le {\eta}_{k}\parallel \nabla {T}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right)\parallel $ |
5: Set [TeX:] $\mathcal {X}^{(k+1)}=\mathcal {X}^{(k)}+\textbf {r}_{k}$ ${\mathcal{X}}^{(k+1)}={\mathcal{X}}^{\left(k\right)}+{\mathbf{r}}_{k}$ |
6: end while |
If the inexact Newton condition (step 4 in Algorithm) is augmented with a sufficiently decreased condition for [TeX:] $\nabla T^{(k)}(\mathcal {X})$ $\nabla {T}^{\left(k\right)}\left(\mathcal{X}\right)$ , the algorithm above can be globally converged:^{38}
Eq. 18
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \Vert \nabla ^{2} T^{(k)}(\mathcal {X}^{(k)}\!+\textbf {r}_{k})\Vert \le [1\!-\!t(1-\eta _{k})]\Vert \nabla T^{(k)}(\mathcal {X}^{(k)}\Vert, t\in (0,1).\nonumber\!\!\!\! \\ \end{eqnarray} \end{document} $$\begin{array}{c}\hfill \parallel {\nabla}^{2}{T}^{\left(k\right)}({\mathcal{X}}^{\left(k\right)}+{\mathbf{r}}_{k})\parallel \le [1-t(1-{\eta}_{k})]\parallel \nabla {T}^{\left(k\right)}({\mathcal{X}}^{\left(k\right)}\parallel ,t\in (0,1).\end{array}$$However, a backtracking strategy is alternatively adopted to transform the augmented condition [Eq. 18] into a practical formulation. At each iteration, an initial inexact Newton step at a specified level is tried, and if it proves unsatisfactory, the inexact Newton steps to higher levels that are solved until a reliable enough step is obtained. The detailed backtracking strategy is shown as:
Table 3
while [TeX:] $\Vert \nabla ^{2} T^{(k)}(\mathcal {X}^{(k)}+\textbf {r}_{k})\Vert \le [1-t(1-\eta _{k})]\Vert \nabla T^{(k)}(\mathcal {X}^{(k)}\Vert$ $\parallel {\nabla}^{2}{T}^{\left(k\right)}({\mathcal{X}}^{\left(k\right)}+{\mathbf{r}}_{k})\parallel \le [1-t(1-{\eta}_{k})]\parallel \nabla {T}^{\left(k\right)}({\mathcal{X}}^{\left(k\right)}\parallel $ do |
Choose θ ∈ [θ_{min}, θ_{max}] |
Update [TeX:] $\textbf {r}_{k}=\theta \textbf {r}_{k}$ ${\mathbf{r}}_{k}=\theta {\mathbf{r}}_{k}$ , η_{k} = 1 − t(1 − η_{k}) |
end while |
Hence, the inexact Newton method with backtracking strategy is incorporated together to enhance convergence from an arbitrary initial guess for unknowns. Up until the end, the proposed reconstruction method is established for the TBI inverse problem. The algorithm shown in Fig. 1 summarizes how the global inexact Newton approach takes advantage of a dynamically sparse regularization technique. The flowchart includes dual-level iterations. The outer iteration controls the whole inexact algorithm, and the sparse regularization term is updated based on the previous iteration during each outer iteration. In order to maintain accuracy of the approximation for [TeX:] $Q_{S}^{(k)}(\mathcal {X}^{(k)})=S(\mathcal {X}^{(k)})$ ${Q}_{S}^{\left(k\right)}\left({\mathcal{X}}^{\left(k\right)}\right)=S\left({\mathcal{X}}^{\left(k\right)}\right)$ , the small ε_{S} is set as 0.02. The inner iterations also fall into two parts: the Krylov solver and the backtracking operation. Here, the preconditioned conjugate gradient is selected as the Krylov solver. p = 1 is taken as an example in Eq. 14.
3.
Results
In this section, an in vivo heterogeneous mouse reconstruction experiment was implemented to demonstrate the feasibility of the proposed method. The experiment was performed on the dual-modality optical/micro-CT in vivo imaging system developed by our group.^{31, 40, 41} The optical detector was a highly sensitive CCD camera (VersArray, Princeton Instruments, Trenton, New Jersey) coupled with a lens (Nikkor, Nikon, Japan). The bioluminescent source was simulated by a home-made luminescent bead, which had an emission spectrum similar to that of a firefly luciferase-based source. Its dimension was about 1.5 mm in diameter and 2.5-mm long. A nude hairless mouse (Nu/Nu, Laboratory Animal Center, Peking University, China) was adopted in this experiment.
Before optical and x ray data acquisition, the CCD was refrigerated to −110°C to reduce dark current noise, and the mouse was anesthetized and the bead was implanted stereotactically into the interspaces between the left and right lobes of the liver. The optical data was collected first. Bioluminescence images from four views were acquired from the mouse surface, with 60 s integration time for each image, and then four corresponding white mouse images were also obtained. After finishing optical acquisition, the mouse was scanned using the micro-computed tomography (CT) to obtain the surface and anatomical structure.
Data processing followed data acquisition. The mouse structure images were reconstructed by the CT cone-beam reconstruction algorithm.^{42} Then, the data were segmented into a heterogeneous volumetric mesh for image reconstruction, as shown in Fig. 2c. This mesh contains 23,752 tetrahedral elements and 4560 discretized nodes with 1092 nodes on the surface. The torso applied for reconstruction covered over 60% of the volume of the mouse body. Then, the demanded photon density projections on the mouse surface were obtained, and the main procedure is summarized in Fig. 2. First, the two sets of data were spatially registered by the corresponding markers on 2D optical images and on 3D CT volume, which is shown in Fig. 2. The image registration results in Fig. 2d were obtained by incorporating Figs. 2a and 2b. Second, the complete-angle outgoing photon density on the mouse surface [Fig. 2f] was projected from the 2D images on the CCD [Fig. 2e] to the 3D surface [Fig. 2c] according to the co-registered image [Fig. 2d]. The source was easily distinguished in CT images, and the actual position of the source could be confirmed at (25.54, 21.31, 8.52), as shown in Fig. 2a. Additionally, the optical properties for each organ were determined with the inverse adding doubling scheme,^{43} as listed in Table 1.
Table 1
Optical properties for each organ in the mouse.
Experimental mouse | Mouse atlas | |||
---|---|---|---|---|
μa | $\mu ^{\prime }_{s}$ μs′ | μa | $\mu ^{\prime }_{s}$ μs′ | |
Heart | 0.022 | 1.129 | 0.058 | 0.963 |
Lungs | 0.071 | 2.305 | 0.195 | 2.173 |
Liver | 0.128 | 0.646 | 0.345 | 0.678 |
Spleen | – | – | 0.345 | 0.678 |
Muscle | 0.075 | 0.586 | 0.086 | 0.429 |
Bone | 0.032 | 2.178 | 0.060 | 2.495 |
3.1.
Reconstruction Comparison
Based on the in vivo mouse experiment, the proposed method was first evaluated by comparing it with other methods. One of the exact Newton methods (as in Algorithm) was selected to demonstrate the predominance of the global inexact Newton method. A recently developed gradient-free method, which was generalized graph cuts, was also employed to further demonstrate the potential and effectiveness of the proposed method.^{31} This method also has the property of global convergence.^{44} As reported previously,^{16, 17, 30, 31} the reconstruction on a large region often results in troublesome problems when using the exact Newton-type method, and a permissible source region is preferred to reduce the number of unknowns and keep reliable reconstruction results.
Figure 3 shows the comparison of the reconstruction results for these methods. The cross sectional images based on the exact Newton method are shown in Figs. 3c and 3(d). The permissable source region in Fig. 3c is {(x, y, z)|22 ⩽ x ⩽ 28, 17 ⩽ y ⩽ 23, 6 ⩽ z ⩽ 10}, and the result in Fig. 3d is reconstructed on the whole region as depicted in Fig. 2c. It is obvious that the selection of the permissable source region is very necessary to obtain reliable reconstruction results. If the permissable region is extended or even no permissable region is adopted, and more unknowns are introduced in the reconstruction, the reconstructed source distribution may deviate further from the real one further, as shown in Fig. 3d. In contrast to this method, as shown in Fig. 3e, both accurate source localization and distribution can be recovered based on the global inexact Newton method applied on the entire region. The reconstruction image in Fig. 3b is the result based on generalized graph cuts. Both of the proposed and generalized graph cuts methods can localize source distribution well compared with the true one. However, the proposed method can also quantify source density, but this was not the case of the gradient-free method. In the reconstructions based on the Newton method and generalized graph cuts, the best regularization parameters were selected and the parameters were in the range of 10^{−5} to 10^{−3}. For the proposed method, the value was 4 × 10^{−2}.
3.2.
Reliability Studies for Imaging Reconstructions
In this part, reconstructions regularized with different orders of magnitude (10^{−1} to 10^{−12}) parameters were performed for in vivo tomographic bioluminescence imaging. As shown in Fig. 4, reconstructions were hardly affected with the regularization parameters, and the bioluminescent source distribution was accurately reconstructed in all cases. Moreover, the results also offered little difference in the way of quantitative information. Figure 5 further demonstrates the nonsensitive performance of the proposed method with a large range of regularization parameters. During four outer iteration steps, the evolution curve of [TeX:] $||\mathcal {MX}-\Phi ^{m}||/||\mathcal {M}^{t}\Phi ^{m}||$ $\left|\right|\mathcal{MX}-{\Phi}^{m}\left|\right|/\left|\right|{\mathcal{M}}^{t}{\Phi}^{m}\left|\right|$ was almost the same with each other, thus only one curve was given here. It appears that the regularization parameters may not affect the tomographic imaging quality very much.
Second the reliability of the proposed method was further validated by changing the initial guess for unknowns [TeX:] $\mathcal {X}^{(0)}$ ${\mathcal{X}}^{\left(0\right)}$ , which were uniformly distributed in a range from 0 to 200. In Fig. 6, based on the results shown above, the bioluminescent source distribution was recovered credibly, and the reconstructions were hardly affected with initials in all cases. Moreover, the results also represented similarly quantitative information with each other. Similar to the curve in Fig. 5, it is also further demonstrated in Fig. 7 that the proposed method can tolerate different initial unknowns and converge reliably. It is also interesting that during iteration, the evolution values of [TeX:] $||\mathcal {MX}-\Phi ^{m}||/||\mathcal {M}^{t}\Phi ^{m}||$ $\left|\right|\mathcal{MX}-{\Phi}^{m}\left|\right|/\left|\right|{\mathcal{M}}^{t}{\Phi}^{m}\left|\right|$ for large initial values would merge into one curve sooner or later, which is almost the same with each other. It is demonstrated that, as proven mathematically,^{38} the global inexact Newton method with dynamic sparse regularization is globally optimized, and the tomographic imaging quality can be maintained.
Moreover, it is noteworthy that all reconstructions were performed on the whole grid, instead of being performed on small permissible source regions. In the mouse experiments, it is not always reliable or feasible to define such regions effectively, so the proposed method is highly applicable to practical tomographic imaging.
3.3.
Evaluation of the Tolerance for Optical Property Mismatch
Based on the advantage of heterogeneous tissue distribution, the optical property in each organ makes an indispensable role for high quality imaging reconstruction. When constructing heterogeneous optical property distribution, a simple and convenient option is to introduce average parameter values from measured or published ranges for individual data.^{30, 45, 46, 47} Although this approach may not be as accurate as straightforward parameter measurements on individual animals, its limitation could be a trade-off for its simplicity and avoidance of additional operations. Here, the reconstruction performance that results from the mismatch in optical properties was evaluated. Two mouse models were employed here, including the above-mentioned experimental mouse and a mouse atlas (which contains 25,783 tetrahedral elements and 4614 nodes with 486 nodes on the surface), as shown in Fig. 8. The atlas was also constructed by our group, and the optical properties for each tissue are listed in Table 1. More details about the atlas can be found in Liu et al.^{31} Based on the mouse atlas, the case with the two sources were considered for evaluating the tolerance for optical property mismatch. For the two models, the optical property mismatch of ±20% and ±50% in both μ_{a} and [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ was introduced into the tissues. Practically, the optical properties may have a combination of positive and negative mismatch of different magnitudes in the absorption and scattering properties. Consequently, four extreme conditions for all of the tissues were considered here: both μ_{a} and [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ were overestimated; both μ_{a} and [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ were underestimated; μ_{a} was overestimated, [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ was underestimated; μ_{a} was underestimated, and [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ was overestimated, as listed in Table 2.
Table 2
The summary for the effect of optical property mismatch of ±20% and ±50% in both μa and $\mu ^{\prime }_{s}$ μs′ . OP. mismatch denotes the bias from real optical properties.
OP. mismatch | |||||||
---|---|---|---|---|---|---|---|
No. | μa | $\mu ^{\prime }_{s}$ μs′ | Reconstruction center | Relative errors | Offset | Artifact | |
1 | +20% | +20% | (25.12, 21.95, 8.32) | (0.42, 0.64, 0.20) | 0.79 | None | |
2 | −20% | −20% | (24.93, 20.94, 8.75) | (0.61, 0.37, 0.23) | 0.75 | None | |
3 | +20% | −20% | (25.61, 22.02, 8.98) | (0.07, 0.71, 0.46) | 0.85 | None | |
4 | −20% | +20% | (25.03, 21.61, 9.02) | (0.51, 0.30, 0.50) | 0.77 | None | |
Model 1 | 5 | +50% | +50% | (26.03, 22.45, 9.82) | (0.49, 1.14, 1.30) | 1.80 | None |
6 | −50% | −50% | (23.86, 20.43, 7.86) | (2.20, 3.01, 1.00) | 2.01 | Exist | |
7 | +50% | −50% | (24.99, 20.91, 8.92) | (0.65, 0.40, 0.40) | 0.86 | None | |
8 | −50% | +50% | (24.75, 20.83, 8.12) | (0.79, 0.48, 0.40) | 1.01 | None | |
1 | +20% | +20% | (21.21, 33.98, 10.81) | (0.79, 0.02, 0.81) | 1.13 | None | |
(21.00, 34.01, 15.11) | (0.50, 0.00, 0.11) | 0.51 | None | ||||
2 | −20% | −20% | (21.51, 33.69, 11.01) | (0.49, 0.31, 1.01) | 1.16 | None | |
(21.00, 33.91, 15.29) | (0.50, 0.09, 0.29) | 0.59 | None | ||||
3 | +20% | −20% | (21.21, 34.01, 10.69) | (0.79, 0.01, 0.69) | 1.06 | None | |
(21.05, 34.00, 15.01) | (0.45, 0.00, 0.01) | 0.54 | None | ||||
4 | −20% | +20% | (21.39, 33.91, 10.83) | (0.61, 0.09, 0.83) | 1.03 | None | |
Model 2 | (21.04, 33.90, 15.31) | (0.46, 0.10, 0.31) | 0.56 | None | |||
5 | +50% | +50% | (21.51, 33.09, 10.11) | (0.49, 0.91, 0.11) | 1.04 | Exist | |
(20.82, 32.29, 13.91) | (0.68, 1.71, 1.09) | 2.14 | Exist | ||||
6 | −50% | −50% | (21.53, 33.82, 11.21) | (0.47, 0.18, 1.21) | 1.31 | Exist | |
(21.01, 33.89, 15.73) | (0.49, 0.11, 0.73) | 0.89 | Exist | ||||
7 | +50% | −50% | (21.28, 34.12, 10.73) | (0.72, 0.12, 0.73) | 1.03 | None | |
(20.83, 33.79, 15.46) | (0.67, 0.21, 0.46) | 0.84 | None | ||||
8 | −50% | +50% | (21.42, 33.91, 10.61) | (0.58, 0.09, 0.61) | 0.86 | None | |
(20.81, 33.79, 15.64) | (0.69, 0.21, 0.64) | 0.96 | None |
The Monte Carlo (MC) method is accurate for the simulation of photon propagation through biological tissues. In the mouse atlas experiment, a MC-based molecular optical simulation environment was employed here to statistically estimate the bioluminescence signal distribution on the mouse surface.^{48, 49} In the simulation settings, two spherical solid sources with a radius of 1.0 mm were located at (22.00, 34.00, 10.00) and (21.50, 34.00, 15.00). A total of 10^{6} photons for each source were tracked for enhancing simulation accuracy.
The reconstruction results are summarized in Table 2, and typical results are shown in Fig. 9. In general, both the mouse and the mouse atlas, the bioluminescent sources were reliably reconstructed. The reconstruction offset is summarized in Fig. 10. In the first case, the maximum location error was 2 mm when there was a –50% mismatch for both μ_{a} and [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ (No. 6 of model 1); in the second case, the maximal error occurred when a +50% mismatch both in μ_{a} and [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ (the second source in No. 5 of model 2) existed with ∼2 mm offset as well. From Table 2, it is observed that the reconstructed sources were localized less ideally with ±50% errors than with ±20%. It consequently seemed that the effects of the mismatch on the tomography results became larger for increasing optical errors. With the mismatch increasing, an artefact also appeared around the source regions [as the red arrow in Fig. 9b], which may be an inevitable side-effect for large optical property errors.^{50}
Another interesting effect seemed that errors in opposite directions partially cancelled each other out, leading to improved source localization. As plotted in Fig. 10, when all tissues had +50% errors in μ_{a} and −50% errors in [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ or −50% errors in μ_{a} and +50% errors in [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ [No. 5, 6 in Figs. 10a and 10b], the localization errors were much better than those with the same sign, and were similar to those with ±20%. The artifacts also appeared to affect the imaging quality only in the same direction with a 50% error (−50% errors for model 1, and ±50% for model 2). Moreover, the reconstructed source became more diffuse with greater errors, compared with lower optical errors. A typical comparison is shown in Figs. 9a and 9b. In other words, when optical errors with the same sign occurred during reconstruction, the imaging quality would be degraded further.
The imaging quality based on the mouse atlas was better than that of mouse experiment. The results in Fig. 9c for +50% errors in μ_{a} and -50% errors in [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ were very similar to that of the +20% errors in μ_{a} and −20% errors in [TeX:] $\mu ^{\prime }_{s}$ ${\mu}_{s}^{\prime}$ in Fig. 9d. The major reason may be the inevitable discrepancies in the optical properties for the experimental mouse tissues. In practice, it is impossible to obtain accurate optical values.^{47, 50} The source in the experimental mouse was located between the left and right lobes of the liver, thus the space was relatively narrow, and the diffuse approximation-based forward model may not fully deal with this special case, which is another affecting factor. Higher order approximation-based models may be helpful for higher imaging quality in narrow structures.^{51, 52}
For all of the reconstructions based on the two mouse models in this part, the regularization parameter was set as 4 × 10^{−2}, and the unknowns were uniformly initialized as zero. It is worth addressing that, based on the mouse atlas using other values of the regularization parameter and the initial unknown guess, the tomography results can be obtained with similar imaging quality as shown in Table 2.
3.4.
Efficiency Studies for Imaging Reconstructions
The high efficiency of the proposed method was also investigated, compared with the exact Newton method, and the generalized graph cuts (GGC) method. As listed in Table 3, four discrete grids of varying sizes of the experimental mouse were utilized.
Table 3
Efficiency comparisons of imaging reconstructions. The size of the grid means the number of points × the number of elements.
No. | Grid size | Newton method | GGC | Proposed method |
---|---|---|---|---|
1 | 2125 × 7761 | 817.86 s | 39.88 s | 20.04 s |
2 | 3048 × 15890 | 2294.66 s | 68.58 s | 28.80 s |
3 | 3714 × 17218 | 4931.42 s | 146.17 s | 42.46 s |
4 | 4560 × 23752 | 7891.52 s | 345.74 s | 82.12 s |
Based on the different methods, the reconstruction time of the four grids is also summarized in Table 3. The reconstruction cost of the global inexact Newton method was much cheaper than the exact Newton one. The efficiency of the proposed method was about 1 to 2 orders-of-magnitude of the exact method. As the grid dimension increased, efficiency predominance became larger, which is also shown in Fig. 11. The key factor is that in the proposed method, only the approximate rather than the exact solution was computed for the next iteration. Moreover, while the efficiency of GGC was also much higher than the exact Newton method, the performance of the proposed method overcame the GGC, which is also clearly depicted in Fig. 11. The reason lay in the complex process for the pairwise terms in the unstructured grid (not based on square pixels or tube vortexes).^{31} On the unstructured grid, the neighborhood for each node was uncertain, and there were more neighbor nodes around a node than the structured grid (based on square pixels or tube vortexes). This process involved a time-consuming operation. Hence, it is seen that the proposed method was very efficient, and was also promising for the inverse problem of TBI. All of the simulations were done on a Intel Core 2 Duo 1.86 GHz PC with 3 GB RAM.
4.
Discussions and Conclusion
In this study, we proposed an efficient global inexact Newton method with a dynamic sparse regularizer for in vivo TBI. In order to make full use of a priori information of the source sparse-distribution, the sparse regularizer is dynamically approximated by a corresponding weighted quadratic norm in each iteration, rather than being adopted explicitly. The inexact Newton with a backtracking optimization technique has the capability to search for the global optimal solution among multiple local minima of the objective function. The in vivo mouse experiment shows that the proposed method can significantly enhance the reconstruction robustness of the regularization parameter and initial values with no a priori assumptions on the bioluminescent source biodistribution. Applying both the experimental data and the Monte Carlo-based mouse atlas data, the sensitivity evaluation for optical property mismatch indicates that even when there exist up to 50% overestimate and underestimate errors in the heterogeneous mouse models, this method can still maintain adequate reconstruction quality. During the reconstructions, it exhibits high computational efficiency and reliable convergence behavior in TBI reconstructions. To summarize, it is demonstrated that the proposed method bears a strong potential to improve the reconstruction capacity for practical tomographic imaging applications.
Generally, data fusion from different modalities plays an indispensable role for high imaging fidelity,^{53, 54, 55} as specially in tomographic bioluminescence imaging, and it is also necessary to combine the information from optical and micro-CT modalities. The data fusion between bioluminescence imaging and micro-CT is helpful to alleviate the ill-posedness and improve imaging quality. First optical images and CT volume can be spatially co-registered, and the captured 2D bioluminescent signal can be projected onto the 3D mouse surface. Second the heterogeneous anatomical structure by micro-CT is naturally fused through the construction of the optical forward model. The fusion imposes such useful a priori information that accurate reconstruction results can be estimated. Moreover, although the anatomical map has been incorporated in the forward model, the heterogeneous priors can also be employed in the inversion process.^{36}
Furthermore, although the diffusion approximation model for TBI reconstructions is very popular, more accurate forward models to describe photon propagation in biological tissues are also necessary for higher imaging quality. Since the emission spectra peaks of the four main luciferase enzymes (Fluc, BGr68, CBRed, and hRluc) are 612, 543, 615, and 480 nm at temperature 37^{ ○}C, respectively, a significant part of the emission spectra deviates from the high-scattering and low-absorbing window. Moreover, when light propagates through a geometrically small volume, the diffusion assumption also becomes invalid to some extent. As a consequence, more complex forward models to compensate the no-diffusion condition have the potential to improve the reconstruction quality further.^{28, 52} After the linear relationship based on improved models has been established, our proposed method in a generalized (l _{p}) regularization framework can be applied for tomographic imaging reconstructions as well.
In conclusion, a dynamically-sparse regularized global method has been presented. Both in vivo experimental and numerical reconstructions have validated that our proposed method has a high capability in maintaining accurate tomographic imaging quality. As discussed above, future work will focus on studying methods for further improving the tomographic imaging performance validated by in vivo experiments with probe-marked tumor models for further research.
Acknowledgments
This study was funded by the National Basic Research Program of China (973 Program) under Grant No. 2011CB707700, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant No. KGCX2-YW-907, the Hundred Talents Program of the Chinese Academy of Sciences, the National Natural Science Foundation of China under Grant Nos. 81027002 and 81071205, and the Science and Technology Key Project of Beijing Municipal Education Commission under Grant No. KZ200910005005.