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

Hongbo Guo; Jingjing Yu; Xiaowei He; Yuqing Hou; Fang Dong; Shuling Zhang

doi:10.1364/BOE.6.001648

1. 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–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₂-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₂-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₂ regularization [15–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₀-norm is an ideal sparsity regularizer, reconstruction based on L₀-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₁-norm based reconstruction methods and produce better results over L₂-norm based methods in scenarios of recovering small localized targets [18–23]. However, solutions of L₁-norm methods are still less sparse and inefficient when heavy-tail 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₀ 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₁-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₁-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.

Method

2.1 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–35]:

{\begin{cases} \nabla \cdot (D_{x} (r) \nabla Φ_{x} (r)) - μ_{a x} (r) Φ_{x} (r) = - Θ_{s} δ (r - r_{l}) \\ \nabla \cdot (D_{m} (r) \nabla Φ_{m} (r)) - μ_{a m} (r) Φ_{m} (r) = - Φ_{x} (r) η μ_{a f} (r) \end{cases}, r \in Ω

where

Φ_{x, m}

denotes the photon flux density, and subscripts x and m denote the excitation and emission light, respectively.

D_{x, m} = 1 / 3 [μ_{a x, a m} + (1 - g) μ_{s x, s m}]

is the diffusion coefficient in biological tissues, with

μ_{a x, a m}

being the absorption coefficient,

g

the anisotropy parameter, and

μ_{s x, s m}

the scattering coefficient.

η μ_{a f} (r)

denotes the unknown fluorescent yield distribution to be reconstructed.

r_{l}

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

\partial Ω

of the domain

Ω

[34]. By solving the diffusion equation with the finite element method, we can obtain a linear model:

A X = Φ

where

X

is an

N \times 1

vector representing the unknown fluorescent source distribution, is an

M \times 1

vector which contains the boundary measurements, and

A

is an

M \times N

weight matrix, that depends on the geometry and the optical properties of the turbid medium.

2.2 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:

\min_{X} {\frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + λ {‖ X ‖}_{P}^{P}}

where

λ

is the regularization parameter and

{‖ X ‖}_{P}^{P}

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,

\min_{X} {\frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + λ {\sum_{i = 1}^{n} | x_{i} |}^{1 / 2}}

and

{\sum_{i = 1}^{n} | x_{i} |}^{1 / 2} = \sum_{i = 1}^{n} (1 / \sqrt{| x_{i} |}) | x_{i} |

, where

| x_{i} |

denotes the absolute value of

x_{i}

. So Eq. (4) was converted into the following form:

\min_{X} {\frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + λ \sum_{i = 1}^{n} (1 / \sqrt{| x_{i} |}) | x_{i} |}

Due to the non-convexity introduced by L_1/2-norm regularizer, it is difficult to solve Eq. (5) directly and efficiently. However, we can approximately transform Eq. (5) into a weighted L₁ minimization problem and solved it iteratively. Let $w_{i} = \frac{λ}{\sqrt{| x_{i}^{t} |}}$ , we get

X^{t + 1} = \arg \min_{X} {\frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + \sum_{i}^{n} w_{i} | x_{i} |}

where

w_{i} > 0

denotes the weight at index

i

, and

x^{t}

denotes the estimated signal in last iteration. To guarantee the feasibility, we used

w_{i} = \frac{λ}{\sqrt{| x_{i}^{t} |} + κ}

for weights calculation, here

κ

is a small positive real number (e.g.

κ = 1 e - 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₁ minimization problem, which can be solved by efficient L₁ 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 $w_{t}$ to $w_{t + 1}$ at the $(t + 1)$ th iteration. To incorporate changes in the weight and quickly compute the new solution of $X$ , we propose the following homotopy program:

\min_{X} {\frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + \sum_{i}^{n} ((1 - ε) w_{t} + ε w_{t + 1}) | x_{i} |}

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]:

a_{i}^{T} (A X^{*} - Φ)=-((1- ε) w_{t} + ε w_{t + 1}) z_{i} for all i \in Γ

| a_{i}^{T} (A X^{*} - Φ) | <(1- ε) w_{t} + ε w_{t + 1}, for all i \in Γ^{c}

where

Γ

denotes the support of

X^{*}

,

z

denotes the sign sequence of

X^{*}

of

Γ

, and

a_{i}^{}

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

\partial X

and the optimality conditions change as [36]:

A_{Γ}^{T} (A X^{*} - Φ) + δ A_{Γ}^{T} A \partial X = - ((1 - ε) W_{t} + ε W_{t + 1}) z + δ (W_{t} - W_{t + 1}) z,

{‖ \underset{p i}{\underset{︸}{a_{i}^{T} (A X^{*} - Φ)}} + δ \underset{d i}{\underset{︸}{a_{i}^{T} A \partial x}} ‖}_{\infty} \leq \underset{q i}{\underset{︸}{(1 - ε) w_{t} + ε w_{t + 1}}} + δ \underset{s i}{\underset{︸}{(w_{t} + w_{t + 1})}}

where

w_{t}

and

w_{t + 1}

denote

| Γ | \times | Γ |

diagonal matrices with their diagonal entries being the values of

w_{t}

and

w_{t + 1}

on

| Γ |

respectively.

With the change of $ε$ , by the new optimality conditions Eq. (9.a) the update direction of $X$ can be written as:

\partial X = {\begin{cases} (^{A_{Γ}^{T} A_{Γ}^{}) - 1} (W_{t} - W_{t + 1}) z o n Γ \\ 0 ...... o t h e r w i s e . \end{cases} .

As we increase $δ$ , the solution moves in the direction $\partial 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_{i \in Γ^{c}} {(\frac{q_{i} - p_{i}}{- s_{i} + d_{i}}, \frac{- q_{i} - p_{i}}{s_{i} + d_{i}})}_{+}, δ^{-} = \min {(\frac{- X_{i}^{*}}{\partial x_{i}})}_{+}

$δ^{+}$ 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 $ε + δ$ , the $X^{*}$ becomes $X^{*} + δ^{*} \partial X$ , 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 $X^{0} = (1, \dots, 1)$ , which means the solution of the first iteration corresponds to a solution of L₁-norm minimization. The iterative reconstruction algorithm stops when the number of iteration reaches a given value k.

3. 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₁-norm methods were selected empirically, ranging from 1e-9 to 1e-14. The maximum iteration number k = 50.

3.1 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].

Table 1. Optical parameters for the mouse organs

View Table | View all tables in this article

A cylindrical fluorescent target with a 0.8mm radius and 1.6mm height was placed in the liver with center at P₀ = (11.9, 6.4, 16.4mm), as shown in Fig. 1(a). The fluorescent yield of the target was set to be. 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.

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 $Z = 16.4 m m$ . 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.

Download Full Size | PDF

We thoroughly compared six methods in single-target reconstruction, including the proposed IRW-L_1/2, L₁-norm based homotopy algorithm (Homotopy-L₁), the incomplete variables truncated conjugate gradient algorithm (IVTCG) [23], TV-norm based two-step iterative shrinkage algorithm [40], L₂-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₁ 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/2method 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.

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.

Download Full Size | PDF

Table 2. Quantitative results in single-target experiment (Best results are in bold)

View Table | View all tables in this article

3.2 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.

Fig. 3 Illustration for simulation settings. (a) shows the tested target positions, where P1, P2 P3, and P4 correspond to Y = 8.4, 10.4, 12.4 and 14.4 mm, respectively. (b) shows different positions of excited nodes, where Z0 = 16.4mm, Z1 = 12.4mm, Z2 = 20.4mm and Z3 = 24.4mm.

Download Full Size | PDF

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₁ have lower FWHM, which means they produce improved spatial resolution over the other methods.

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.

Download Full Size | PDF

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₁ takes the second place in most of the test cases. Therefore, only IRW-L_1/2 and Homotopy-L₁ are considered in the following simulations.

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.

Download Full Size | PDF

3.3 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₁ and IRW-L_1/2 algorithm suffer from some degree of fluctuation, the performance of IRW-L_1/2 is more stable than L₁-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₁.

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.

Download Full Size | PDF

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₁-norm method. However, IRW-L_1/2 is less sensitive to reduction of measurements and outperforms Homotopy-L₁. 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₁ deteriorates rapidly, with location deviation more than 1mm.

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.

Download Full Size | PDF

The stability analysis in this section demonstrates that the performance of the proposed reconstruction algorithm IRW-L_1/2 is superior to Homotopy-L₁. 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.

3.4 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₁ 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₁, the related result has larger deviation in target positions and fluorescence yield.

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₁; (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.

Download Full Size | PDF

Table 3. Quantitative results for reconstruction of two targets with 4.5mm separation

View Table | View all tables in this article

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(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₁. Table 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₁ fail to discriminate such close two targets. It means that IRW-L_1/2 provides better spatial discrimination than Homotopy-L₁.

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₁. (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.

Download Full Size | PDF

Table 4. Quantitative results for reconstruction of two targets with 3mm separation

View Table | View all tables in this article

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.

3.5 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 [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.

Fig. 10 Reconstructed results of in vivo experiment on BALB/C mouse with implanted fluorophore. (a) shows the forward mesh and the measurements mapped on surface, (e) shows the segmented torso with implanted fluorophore in a glass tube; (b) and (f) are the results of IRW-L_1/2; (c) and (g) are the results of Homotopy-L₁; (d) and (e) are the results of IVTCG. Subfigures (b), (c) and (d) correspond to the cross-section views, where the small black circle denotes the true source; (f), (g) and (h) are the corresponding 3D view, where the red cylinder is the true target while the green zone denotes the reconstructed one.

Download Full Size | PDF

Table 5. Optical parameters of the mouse organs at 670 and 710 nm

View Table | View all tables in this article

The IVTCG algorithm was used in [23] to reconstruct the implanted fluorophore, which is also an L₁-norm based reconstruction method. Table 6 presents detailed reconstruction results by IRW-L_1/2, Homtopy-L₁ and IVTCG. Figure 10 shows the experimental setting, forward photon distribution on surface and the corresponding reconstruction results in cross-section views and 3D views, where the color bar of RFY is normalized. To illustrate the reconstruction accuracy, we plot the transverse views of the IRW-L_1/2 result in Fig. 11 and compared them with the corresponding CT slices. As we can see in Fig. 10 and Table 6, two L₁-norm algorithms (Homtopy-L₁ and IVTCG) produce very similar results with a deviation of 1.73mm, while the proposed IRW-L_1/2 algorithm obtain better reconstruction with a deviation of 1.30mm. Due to the PNZ of IRW-L_1/2 is 0.22%, which is smaller that of L₁-norm methods, the spatial distribution in Fig. 10(f) looks closer to the true target than that in Fig. 10(g) and Fig. 10(h). Detailed comparisons are summarized in Table 6. In general, compared to Homtopy-L₁ and IVTCG, our IRW-L_1/2 algorithm produced more confined solution with smaller LE, clean background and higher RFY.

Table 6. Quantitative results of in vivo experiment

View Table | View all tables in this article

Fig. 11 Transverse views of the result by IRW-L_1/2 and the comparison with the corresponding CT slices. The area surrounded by the green lines denotes the true target determined by CT.

Download Full Size | PDF

4. Discussion and conclusion

To accurately recover small fluorescent targets, we combined L_1/2-norm regularization in FMT reconstruction. By iterative reweighting, we transformed the nonconvex L_1/2-norm problem into a series of weighted L₁-norm minimization problems and efficiently solved them with a homotopy-based algorithm IRW-L_1/2. We demonstrated the performance of IRW-L_1/2 in FMT with simulated data on a heterogeneous mouse model and in vivo experimental data acquired by a dual-modality FMT/Micro-CT system. By thorough comparison study on single-target reconstruction, we show the proposed IRW-L_1/2 algorithm outperformed several typical reconstruction algorithm based on L₁-norm, L₂-norm, TV-norm and ART. The reconstruction results of IRW-L_1/2 also revealed the sparsity-enhancing ability of L_1/2-norm regularization, which is proved to be useful in resolving multiple targets and improving reconstruction quality and spatial resolution. Double-target simulations results in section 3.4 show that IRW-L_1/2 successfully recovered two targets with 3mm separation, while Homtopy-L₁ algorithm produced inferior results in this challenging situation.

The stability of IRW-L_1/2 over noise and reduction of excitation nodes was also analyzed with simulations. When the noise level in measurement gradually increased, we observed the location accuracy of the proposed algorithm kept steadily and the reconstructed fluorescent yield slightly fluctuated. In addition, although the reduction of and rec excitation nodes aggravated the ill-posedness of FMT reconstruction, the performance of IRW-L_1/2 was still satisfactory. Even in the extreme case where only one excited node was used, the proposed method still successfully recovered the target with a deviation of 0.81mm, which exhibited potential of IRW-L_1/2 in FMT reconstruction with one view measurements.

Although the reconstructed fluorescent yield of our algorithm was generally larger than that of the comparative algorithms, it was still smaller than the true value. New strategies that can further improve fluorescent yield will be our future concern.

In conclusion, both numerical simulations and in vivo experiment validated the feasibility and efficiency of L_1/2-norm regularization and IRW-L_1/2 algorithm in FMT reconstruction. Comparative studies between IRW-L_1/2 algorithm and L₁-norm algorithms demonstrated that IRW-L_1/2 is an alternative method to improve FMT reconstruction quality.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant Nos. 61372046and 61401264), the Research Fund for the Doctoral Program of Higher Education of China (New Teachers) (Grant No. 20116101120018), the China Postdoctoral Science Foundation Funded Project (Grant No. 2012T50814), Science and Technology Plan Program in Shaanxi Province of China(Grant No. 2012 KJXX-29).

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

3. X. Montet, V. Ntziachristos, J. Grimm, and R. Weissleder, “Tomographic fluorescence mapping of tumor targets,” Cancer Res. 65(14), 6330–6336 (2005). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

7. K. Radrich, P. Mohajerani, J. Bussemer, M. Schwaiger, A. J. Beer, and V. Ntziachristos, “Limited-projection-angle hybrid fluorescence molecular tomography of multiple molecules,” J. Biomed. Opt. 19(4), 046016 (2014). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef]

29. D. Zhu and C. Li, “Nonconvex regularizations in fluorescence molecular tomography for sparsity enhancement,” Phys. Med. Biol. 59(12), 2901–2912 (2014). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef] [PubMed]

33. A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

44. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [CrossRef] [PubMed]

Material	$μ_{a x} ({mm}^{- 1})$	$μ_{s x}^{'} ({mm}^{- 1})$	$μ_{a m} ({mm}^{- 1})$	$μ_{s m}^{'} ({mm}^{- 1})$
Muscle	0.0052	1.08	0.0068	1.03
Heart	0.0083	1.01	0.0104	0.99
Lungs	0.0133	1.97	0.0203	1.95
Liver	0.0329	0.70	0.0176	0.65.
Kidneys	0.0660	2.25	0.0380	2.20
Stomach	0.0114	1.74	0.0070	1.36

Method	True Center (mm)	Reconstructed Center (mm)	LE (mm)	RFY (mm−1)	QE (%)	PNZ (%)
IRW-L1/2	(11.9,6.4,16.4)	(11.6,6.5,16.8)	0.54	0.038	24	1.21
Homotopy-L1	(11.9,6.4,16.4)	(11.9,6.8,15.7)	0.81	0.026	48	1.32
IVTCG	(11.9,6.4,16.4)	(12.1,6.7,14.8)	1.67	0.004	89	3.25
TV	(11.9,6.4,16.4)	(10.6,6.4,15.0)	1.92	0.004	89	4.01
L2	(11.9,6.4,16.4)	(13.2,5.0,15.6)	2.11	0.004	89	8.54
ART	(11.9,6.4,16.4)	(10.4,5.8,16.1)	1.69	0.001	98	14.23

Method	Fluorescent Target	Actual Center (mm)	Reconstructed Center (mm)	LE (mm)	RFY (mm⁻¹)	QE (%)	PNZ (%)
IRW-L_1/2	T1	(11.8, 6.4,16.4)	(11.3, 6.3,15.8)	0.78	0.031	38	1.59
IRW-L_1/2	T2	(11.8,10.9,16.4)	(11.7,10.1,15.8)	1.02	0.032	36	1.59
Homotopy-L₁	T1	(11.8,6.4,16.4)	(12.5, 7.6,16.2)	1.56	0.021	59	1.72
Homotopy-L₁	T2	(11.8,10.9,16.4)	(11.0,13.2,16.0)	2.54	0.010	80	1.72

Method	Fluorescent Target	Actual Center (mm)	Reconstructed Center (mm)	LE (mm)	RFY (mm⁻¹)	QE (%)	PNZ (%)
IRW-L_1/2	T1	(11.8,10.8,16.4)	(11.4,11.5,15.6)	1.12	0.033	34	1.23
IRW-L_1/2	T2	(11.8,13.8,16.4)	(11.8,14.6,15.6)	1.12	0.044	3	1.23
Homotopy-L₁	T1	(11.8,10.8,16.4)	(11.8,12.8,15.6)	2.32	0.022	56	2.24
Homotopy-L₁	T2	(11.8,13.8,16.4)	(12.6,13.7,17.0)	0.99	0.046	0.8	2.24

Material	670mm		710mm
Material	$μ_{a x} ({mm}^{- 1})$	$μ_{s x}^{'} ({mm}^{- 1})$	$μ_{a m} ({mm}^{- 1})$	$μ_{s m}^{'} ({mm}^{- 1})$
Muscle	0.075	0.412	0.043	0.350
Heart	0.051	0.944	0.030	0.870
Lungs	0.170	2.157	0.097	2.093
Liver	0.304	0.668	0.176	0.626
Kidneys	0.058	2.204	0.034	2.021

Improved sparse reconstruction for fluorescence molecular tomography with L_1/2 regularization

Abstract

1. Introduction

Method

2.1 Photon propagation model

2.2 The framework of L_1/2-norm regularization

3. Experiments and results

3.1 Single-target reconstruction on 3D digital mouse atlas

3.2 Influence of depth and relative position of the excitation nodes on algorithm performance

3.3 Stability analysis

3.4 Reconstruction of double-target

3.5 In vivo evaluation with implanted fluorophore

4. Discussion and conclusion

Acknowledgments

References and links

Cited By

Figures (11)

Tables (6)

Equations (13)

Biomedical Optics Express

Method	True Center (mm)	Reconstructed Center (mm)	LE (mm)	RFY (mm^-1)	QE (%)	PNZ (%)
IRW-L_1/2	(28.5,25.7,13.4)	(29.6,25.3,13.9)	1.30	0.021	47.5%	0.22
Homotopy-L₁	(28.5,25.7,13.4)	(28.3,27.4,12.9)	1.73	0.016	59.7%	0.55
IVTCG	(28.5,25.7,13.4)	(28.3,27.4,12.9)	1.73	0.014	64.9%	0.69

Abstract

1. Introduction

Method

2.1 Photon propagation model

2.2 The framework of L1/2-norm regularization

3. Experiments and results

3.1 Single-target reconstruction on 3D digital mouse atlas

3.2 Influence of depth and relative position of the excitation nodes on algorithm performance

3.3 Stability analysis

3.4 Reconstruction of double-target

3.5 In vivo evaluation with implanted fluorophore

4. Discussion and conclusion

Acknowledgments

References and links

Cited By

Figures (11)

Tables (6)

Equations (13)

Biomedical Optics Express

2.2 The framework of L_1/2-norm regularization