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,¹⁶ which is given by

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:

Based on the finite element theory, Eqs. 1, 2 are discretized by piecewise linear finite elements, and a matrix-vector equation is integrally assembled.³³ 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:

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:

Now let us begin to construct the dynamic sparse regularized function. Instead of solving the sparse l _p problem,²⁰ a quadratic version is generated to approximate the sparse term for each iteration.³² In order to replace the l _p term by the l ₂ one and dynamically regularize the objective function, the quadratic function is defined as:

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,³⁷ [TeX:] $\tau _{S,\epsilon _{S}}$ ${\tau }_{S,{\epsilon }_S}$ (for some small ε_S) is defined as

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:

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:

The corresponding gradient and Hessian matrix are easily derived

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:

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:³⁸

Table 1

Newton method.

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.³⁹ 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.

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:³⁸

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

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.

Fig. 1

The algorithmic flowchart for the proposed imaging reconstruction method.

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.³¹ This method also has the property of global convergence.⁴⁴ 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⁻⁵ to 10⁻³. For the proposed method, the value was 4 × 10⁻².

Fig. 3

The comparison of cross sectional images for reconstructed bioluminescence source distribution. (b) The result of generalized graph cuts method on the whole reconstruction region. (c) and (d) The result of the Newton method on a small permissable source region and the whole reconstruction region, respectively. (e) The result of the proposed method on the entire reconstruction region. (a) and (f) are the corresponding CT slices. Note that the colorbar scales on the right vary with each reconstruction.

3.2.

Reliability Studies for Imaging Reconstructions

In this part, reconstructions regularized with different orders of magnitude (10⁻¹ to 10⁻¹²) 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.

Fig. 4

Reconstruction results with different regularization parameters. The center is the corresponding CT slice. All of the unknowns are set to be 0 for all cases.

Fig. 5

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|$ as a function of the iteration steps with initial unknowns of 0 uniformly. Four outer iteration steps are undergone during reconstruction. Note that the curve with different regulation parameters is very similar, hence just one curve is plotted here.

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,³⁸ the global inexact Newton method with dynamic sparse regularization is globally optimized, and the tomographic imaging quality can be maintained.

Fig. 6

Reconstruction results with different initial guesses. The initials are uniform for all unknowns. The center is the corresponding CT slice. The regularization parameter in all cases is selected as 4 × 10⁻².

Fig. 7

The evolution curve as a function of iteration steps with different initial unknowns. Four outer iteration steps are undergone.

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.³¹ 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.

Fig. 8

The volumetric mesh of the heterogeneous mouse atlas, including heart, lungs, liver, spleen, muscle, and bone.

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.

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⁶ 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.⁵⁰

Fig. 9

The reconstruction results for both models when an optical property mismatch is introduced. (a) and (b) are results based on the experimental mouse, and (c) and (d) are based on the mouse atlas. (a) and (c) The results with +20% mismatch of both μ_a and [TeX:] $\mu ^{\prime }_{s}$ ${\mu }_s^{\prime }$ . (b) and (d) The results with mismatch of +50% μ_a and −50% [TeX:] $\mu ^{\prime }_{s}$ ${\mu }_s^{\prime }$ . In (b), the red arrow points to the artifact around the real source.

Fig. 10

Comparison for the source localization offset with different optical property mismatch. (a) The localization offset of the experimental mouse. (b) The localization offset for two sources based on the mouse atlas.

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⁻², 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.

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).³¹ 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.