An adaptive support driven reweighted L1-regularization algorithm for fluorescence molecular tomography

Junwei Shi; Fei Liu; Huangsheng Pu; Simin Zuo; Jianwen Luo; Jing Bai

doi:10.1364/BOE.5.004039

1. Introduction

Fluorescence molecular tomography (FMT) is a non-invasive biomedical imaging modality that employs near infrared light to monitor the three dimensional (3D) cell and molecular activities inside biological tissues [1,2]. In the past decade, various studies have been carried out with FMT at the molecular level, including cancer imaging, drug delivery, stem cell trafficking, enzyme activity monitoring, etc [3–6].

However, due to the high scattering of photons in tissues and the relatively small amount of available measurements, the inverse problem of FMT reconstruction is highly ill-posed and ill-conditioned. Moreover, the reconstruction of FMT is highly susceptible to the noise and model errors and faces various challenges on quantitation and spatial resolution. In order to obtain high-quality FMT reconstruction results, different regularization methods have been proposed to alleviate the effect of ill-posedness. Among these methods, Tikhonov as a famous L2-regularization method is widely used [6]. Although L2 regularization is simple to implement and can be efficiently solved, it promotes smooth solutions and thereby suppresses the high-frequency components, resulting in blurred results with low spatial resolution. In most in vivo applications, the fluorophore distribution is usually sparse because the fluorescent probes are designed to accumulate in relatively small regions such as in tumors. Hence, L1-regularization methods are extensively studied and found to yield reconstruction results with high spatial resolution [7–9]. And total variation (TV) regularization as one kind of sparsity regularization was also employed in FMT to promote local smoothness [10]. Furthermore, the method combining L1 and TV regularization was also brought into FMT reconstruction to preserve the boundary and promote sparsity [11].

Both L2 and L1 regularization methods in the minimization scheme result in convex optimization. Recently, in many fields (such as FMT [12,13], compressed sensing [14–17] and diffuse optical tomography (DOT) [18]), non-convex regularization has attracted lots of attention and is found to have greater power than L1 regularization in enhancing sparsity. The non-convex regularization is usually solved using the iteratively reweighted L1-minimization method, where the weights are estimated as the reciprocal value of the previous solution. In order to avoid dividing zero-value component, a positive number as a stability parameter is needed in the denominator [14]. However, the choice of this stability parameter affects the stability of the solution. Based on the support estimate, the support driven reweighted L1-minimization method can avoid the use of aforementioned stability parameter [15,16].

In this paper, an adaptive support driven reweighted L1-regularization (ASDR-L1) algorithm is proposed to solve the non-convex regularization problem for FMT. ASDR-L1 is related to the restarted L1 regularization-based nonlinear conjugate gradient (re-L1-NCG) algorithm previously proposed in Refs [7,8]. In ASDR-L1, the original L1-minimization problem is replaced by a sequence of reweighted L1-minimization sub-problems with spatially updated weights applied to the adaptive support estimate. The ASDR-L1 algorithm has two integral parts: an adaptive support estimate and the iteratively updated weights. Through the proposed threshold strategy, the support estimate can be adaptively obtained in each outer iteration of ASDR-L1. The updated weights for the next iteration are a function of the preceding iterative solution in the corresponding support estimate. It should be noted that the spatially dependent weights in ASDR-L1 means that there are different regularization parameters at different spatial locations. Thus ASDR-L1 can be considered as a kind of spatially variant regularization algorithms as proposed in Ref [19]. In order to evaluate the performance of the proposed algorithm, physical phantom and in vivo mouse experiments were carried out. The results demonstrate that the proposed ASDR-L1 algorithm outperform the Tikhonov and re-L1-NCG methods in terms of relative quantitation and spatial resolution, when the fluorophores are located at different depths.

The structure of the paper is organized as follows. In section 2, the reconstruction methodology is presented. Firstly, the mathematical framework of the forward model for FMT is briefly described. Secondly, the ASDR-L1 algorithm is formulated. Finally, the quantitative metrics is proposed to evaluate the quality of the results reconstructed by different methods. In section 3, physical phantom and in vivo mouse experiments are performed to illustrate the performance of the proposed algorithm. In section 4, we discuss relative issues and conclude our work.

2. Methods

2.1 Forward model for FMT

FMT aims to discover the 3D distribution of the 〉uorescence source inside the object through the forward model which predicts light propagation in biological media. The photon propagation can be modeled using the integro-differential equation known as the radiative transfer equation (RTE). Owing to the prodigious computational burden of RTE, the forward model used to predict photon propagation in highly scattering media is usually approximated as the coupled diffusion equations (DE) with the Robin-type boundary condition [20]. In this paper, the normalized Born approximation to DE is adopted to linearize the nonlinear FMT problem [21]. This normalization scheme can offer significant experimental advantages since it is independent of incident light intensity, detector gains and correct instrumentation response function. Furthermore, it is insensitivity to the theoretical errors including modeling of the boundary conditions and knowledge of the exact tissue optical properties. The normalized Born average intensity $ϕ_{f} (r_{d}, r_{s})$ at location $r_{d}$ corresponding to an illumination spot located at $r_{s}$ can be formulated as follows [21]:

ϕ_{f} (r_{d}, r_{s}) = \frac{ϕ_{e m} (r_{d}, r_{s})}{ϕ_{e x} (r_{d}, r_{s})}

where

ϕ_{e m} (r_{d}, r_{s})

and

ϕ_{e x} (r_{d}, r_{s})

denote the average intensity at emission wavelength

λ_{2}

and excitation wavelength

λ_{1}

, respectively. Correspondingly, the normalized Born approximation can be formulated as follows [21]:

ϕ_{f} (r_{d}, r_{s}) = \frac{S_{0} υ}{G (r_{s}, r_{d}, λ_{1}) D^{f}} \int G (r_{d}, r, λ_{2}) x (r) G (r_{s}, r, λ_{1}) d^{3} r

where

G (r_{s}, r, λ_{1})

denotes the analytically calculated photon field at excitation wavelength

λ_{1}

induced at position

r

by a source at position

r_{s}

in the tissue.

G (r_{d}, r, λ_{2})

is the Green’s function which describes photon propagation from point

r

to the detector point

r_{d}

at emission wavelength

λ_{2}

.

D^{f}

is the diffusion coefficient of tissue at the emission wavelength

λ_{2}

,

υ

is the speed of light in the tissue, and

S_{0}

is a calibration factor accounting for system gain and attenuation factors. Kirchhoff approximation (KA) as an analytical method is used to solve the FMT forward model, in order to ease the computation burden [22]. After the image domain is discretized, the FMT problem can be formulated as the following linear matrix equation:

I^{k}

where W denotes the sensitive matrix mapping the unknown vector of fluorophore concentrations X into the vector of normalized Born average intensity

Φ_{f}

. Detailed descriptions can be found in Refs [20,21].

2.2 FMT reconstruction based on reweighted L1 regularization

FMT recovers the internal 〉uorophore distribution by minimizing the discrepancy between measurements from the FMT system and predicted values from the forward model. No direct solutions for Eq. (3) exist, as the sensitive matrix W is severely ill-conditioned. In the practical applications, the fluorophore distribution is sparse in the entire reconstruction domain. Thus sparsity regularization as one kind of a priori information can be used to guarantee the uniqueness and stability of the FMT inverse problem [9]. Based on the sparsity regularization for FMT, the corresponding optimization function is formulated as the following L0-norm regularization form:

\min_{X \geq 0} {‖ W X - Φ_{f} ‖}_{2}^{2} + λ {‖ X ‖}_{0}

where L0 norm

{‖ X ‖}_{0}

denotes the number of nonzero values in X. Although Eq. (4) is the direct formulation for the sparsity regularization, this optimization problem is generally impossible to be solved as its solution usually requires an intractable combinatorial search. Because L1 norm is the closest convex approximation to L0 norm, L1 regularization is usually employed to obtain the suboptimal solution of Eq. (4). However, as a kind of non-convex approximation to L0-norm regularization, the reweighted L1 regularization was reported to outperform the regular L1 regularization [14].

Based on the reweighted L1 regularization, the FMT reconstruction modeled as a regularized least squares problem with a non-negativity constraint is formulated as follows:

\min_{X \geq 0} {‖ W X - Φ_{f} ‖}_{2}^{2} + λ {‖ R^{k} X ‖}_{1}

where

λ

is a regularization parameter. And the positive diagonal matrix R^k, which is updated in the k-th iteration, denotes the weighted matrix. By comparing Eq. (4) with Eq. (5), it can be seen that in the reweighted L1 regularization the iteratively updated L1 norm

{‖ R^{k} X ‖}_{1}

is used to approximate the L0 norm

X_{L 1 - N C G}^{k}

. Thus, unlike the convex L1 regularization, the overall reweighted L1 regularization is non-convex. However, in order to obtain the solution of the reweighted L1-regularization problem, the L1-regularizaed sub-problem is needed to be reweighted and solved in the k-th iteration.

In order to roughly level off the differences in detection sensitivities, the columns of the sensitive matrix W defined in Eq. (3) are normalized [7,8]. Correspondingly, the resultant preconditioned optimization function is rewritten as follows:

X_{r e - L 1} = \arg \min_{X \geq 0} (f (X) = {‖ W_{n e w} X - Φ_{f n e w} ‖}_{2}^{2} + λ {‖ R^{k} X ‖}_{1})

W_{n o r} (i, j) = {\begin{matrix} {‖ W_{i} ‖}_{2}^{- 1} i = j \\ 0 \begin{matrix}  \end{matrix} i \neq j \end{matrix}

W_{n e w} = W W_{n o r}; Φ_{f n e w} = Φ_{f} / \max (Φ_{f})

where W_i is the i-th column of matrix W. In this work, the ASDR-L1 algorithm is proposed to solve Eq. (6). ASDR-L1 is a modification of the previously proposed re-L1-NCG algorithm in Refs [7,8]. Like re-L1-NCG, in ASDR-L1, the L1-regularized nonlinear conjugate gradient (L1-NCG) method with backtracking line search is used to solve the L1-regularized sub-problem [7]. The main steps of the proposed algorithm are summarized below.

Algorithm 1 Adaptive Support Driven Reweighted L1-regularization (ASDR-L1)

Parameters setting: Set

k_{i n n e r} = 100

as the maximum inner iterations for L1-NCG,

k_{o u t e r}

as the maximum outer iterations for ASDR-L1. k = 1. The other parameters of L1-NCG can be found in Ref [7]. Un-weighted L1-minimization: For Eq. (6) without non-negative constraint, compute the solution

X_{L 1 - N C G}^{1}

using L1-NCG with

k_{i n n e r}

iterations, where the weight matrix R¹ in Eq. (6) is an identity matrix. Reweighted L1-minimization iteration (RLI): (1) In order to achieve the support estimate for ASDR-L1, set

x_{i, L 1 - N C G}^{k} = 0

when

| x_{i, L 1 - N C G}^{k} | < 0.01 \times \max (X_{L 1 - N C G}^{k})

.

x_{i, L 1 - N C G}^{k}

denotes the i-th element of

X_{L 1 - N C G}^{k}

. (2) Achieve the support estimate I^k = {indices with the nonzero values of

X_{L 1 - N C G}^{k}

}.

{\tilde{X}}_{}^{t e m} = n o n z e r o (X_{L 1 - N C G}^{k})

;

{\tilde{X}}_{}^{R} = \max ({\tilde{X}}_{}^{t e m}, 0)

;

M a x = \max ({\tilde{X}}_{}^{R}) = \max ({\tilde{X}}_{}^{t e m})

(3) Generate the weighted matrix R^k for the reweighted L1 regularization.

R^{k} = \max ({\tilde{X}}_{}^{t e m}) \times {[d i a g (g ({\tilde{X}}_{}^{t e m}, M a x)]}^{- 1} = [\begin{matrix} \frac{M a x}{g ({\tilde{x}}_{1}^{t e m}, M a x)} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & \frac{M a x}{g ({\tilde{x}}_{n}^{t e m}, M a x)} \end{matrix}]

g ({\tilde{x}}_{i}^{t e m}, M a x) = {\begin{matrix} {\tilde{x}}_{i}^{t e m} \begin{matrix}  \end{matrix} i f \begin{matrix}  \end{matrix} {\tilde{x}}_{i}^{t e m} > 0 \\ M a x \begin{matrix}  \end{matrix} i f \begin{matrix}  \end{matrix} {\tilde{x}}_{i}^{t e m} < 0 \end{matrix}

(4) In the matrix

W_{n e w}

, remove columns corresponding to zero values of

X_{L 1 - r e}^{k}

and obtain the new matrix

W_{n e w}^{r e}

. Then the optimization function turns into:

X_{L 1 - N C G}^{k + 1} = \arg \min_{X \geq 0} (f (X) = {‖ W_{n e w}^{r e} X - Φ_{f n e w} ‖}_{2}^{2} + λ {‖ R^{k} X ‖}_{1})

(5) Solve the new optimization problem for

X_{L 1 - N C G}^{k + 1}

using L1-NCG with

k_{i n n e r}

iterations, where the initial value is set as

{\tilde{X}}_{}^{R}

and the search direction is restarted as the negative gradient of

f (X)

.(6)

| X_{L 1 - N C G}^{k} |

. If

k > k_{o u t e r}

or

s i z e (I^{k}) - s i z e (I^{k - 1}) < 0.3 * s i z e (I^{k})

, quit the iteration.Output:

X_{r e - L 1}

denotes the reweighted L1-regularized solution to Eq. (6). In

X_{r e - L 1}

, the values corresponding to the indices in

I^{k}

is

X_{L 1 - N C G}^{k}

, and the rest values are zero.

The solution vector $X_{L 1 - N C G}^{1}$ obtained from the un-weighted L1-minimization is used to construct the first support estimate and the initial values for the following loop called as reweighted L1-minimizaton iteration (RLI).

In the first two steps of RLI, support estimate is generated by thresholding the absolute solution vector $| X_{L 1 - N C G}^{k} |$ , in order to avoid deleting the correct support set existing in the indices corresponding to the negative values. This threshold strategy aims to avoid dividing zero-value component in the solution vector, when constructing the weighted matrix for the reweighted L1 regularization. In general, the recover process tends to be reasonably robust to the choice of this threshold which is set to 0.01 in this work.

In the third step of RLI, the weighted matrix R^k for the reweighted L1 regularization is constructed based on the adaptive support estimate achieved in the first two steps of RLI. Based on the formula for R^k, the updated regularization parameters at different locations can be realized. As mentioned in Ref [23], a large regularization parameter in the weighted L1-minimization sub-problem makes the results concentrate on a small number of large values, whereas a small regularization parameter tends to make the values evenly distributed. As the value of regularization parameter has effect on the resolution and contrast in the reconstruction result [23], in order to realize depth compensation for FMT [23], larger reconstruction results will be assigned with smaller regularization parameters and vice versa, which is the basis to construct matrix R^k. In matrix R^k, the weights are normalized by the maximum in the support estimate, in order to realize relative quantitation. And in view of the non-negative constraint of Eq. (6), different weighted strategies are assigned to the positive and negative solutions, respectively. When the solution in the support estimate is positive, the corresponding weight is set to $M a x / {\tilde{x}}_{i}^{t e m}$ , whereas the weight corresponding to the negative solution is set to 1 = Max/Max. On one hand, this weighted strategy for the negative solutions will not magnify the influence of negative solutions on the next iteration. On the other hand, this unified weight (1 = Max/Max) can help to pick up the support existing in the indices corresponding to the negative solutions.

In the fourth step of RLI, by means of the constructed weighted matrix R^k, the corresponding L1-regularized sub-problem is reweighted in the support estimate. And in the fifth step, the reweighted L1-minimization sub-problem is solved using L1-NCG with the updated initial value and the restarted search direction.

The sixth step of RLI is used to determine whether or not it is time to terminate the proposed method. In this work, the proposed ASDR-L1 algorithm is terminated when the relative size reduction between two consecutive support estimates is below a certain threshold, because the reconstruction result will be too sparse if this size reduction is too small. In addition, the proposed algorithm will be terminated when the maximum outer iteration number k_outer is reached. In this paper, the parameter k_outer is manually optimized as 10.

In this work, the Tikhonov method as a widely used analytical reconstruction method is adopted for comparison, where the suboptimal regularization parameter is experientially selected based on $10^{- 2} \times t r (W_{n e w} W_{n e w}^{T})$ . In addition, re-L1-NCG as an L1-regularized method is adopted to solve Eq. (6), where R^k is set as the identity matrix. Because ASDR-L1 and re-L1-NCG have the similar scheme, the comparison between these two methods can well demonstrate the performance differences between reweighted L1 regularization and L1 regularization, where the regularization parameters for these two methods are set to the same value of 15 (as described in [7,8]).

2.3 Quantitative metrics for FMT reconstruction

In order to quantitatively analyze the performance of the algorithms, the relative quantitation index RQ is employed [24]:

R Q = \frac{| x_{1}^{\max} - x_{2}^{\max} |}{| x_{1}^{\max} + x_{2}^{\max} | / 2}

where

x_{1}^{\max}

and

x_{2}^{\max}

are the maximum reconstructed intensities for each target. And these values come from the profile along a given line crossing the reconstruction results on the cross-sectional images. Because absolute quantitation for FMT depends on many factors which cannot be accurately measured in practice, relative quantitation is adopted here. When two identical targets are embedded in the tissue, the index RQ denotes the reconstruction performance in term of relative quantitation. And a smaller RQ indicates better relative quantitation performance.

3. Experiments and results

3.1 Physical phantom experiments

In order to evaluate the performance of the proposed algorithm in terms of relative quantitation and spatial resolution, physical phantom experiments were conducted based on the experimental system previously developed in our laboratory [25]. The schematic of this system is illustrated in Fig. 1. In the experiments, two identical fluorophores with different edge-to-edge distances (EED) were embedded at different depths in the phantom. The phantom was made of a transparent glass cylinder with a diameter of 3 cm and height of 6 cm. The cylinder filled with 1% intralipid (with absorption coefficient of 0.02 cm⁻¹ and reduced scattering coefficient of 10 cm⁻¹) was employed as the phantom, where the optical properties are nearly homogeneous. The fluorescence targets in the phantom were two cylinders (with a diameter of 0.4 cm) containing 20 μL indocyanine green (ICG) with a concentration of 1.3 μM. In order to excite ICG, a band-pass filter with a center wavelength of 770 nm and full width half maximum (FWHM) of 10 nm was used in front of a Xenon lamp (Asahi Spectra, Torrance, CA, USA). Through the focusing lens, point incident light was generated to excite the phantom. In order to collect the fluorescence light, a band-pass filter with a center wavelength of 840 nm and FWHM of 12 nm was employed in front of the high sensitive cooled charge coupled device (CCD) camera (iXon DU-897, Andor Technologies, Belfast, Northern Ireland). In order to collect the excitation light, a neutral density filter with an optical density of 2 was used in front of the CCD camera. The detector field of view (FOV) corresponding to one point excitation source was about 130° and 36 projections were adopted in both fluorescence and excitation light acquisition.

Fig. 1 Sketch of the free-space full-angle FMT system.

Download Full Size | PDF

Through discretizing the phantom, Eq. (3) was obtained. The KA method was used to solve the FMT forward model. The size of the sensitivity matrix W was determined by the product of the number of measurements utilized and the number of voxels employed to discretize the volume of interest. In the physical phantom experiments, the size of the sensitivity matrix W was 9567 × 15133, where the detector spatial sampling was set as 0.25 cm and the discretized mesh resolution of the phantom was set as 0.11 × 0.11 × 0.11 cm³. The normalized Born measurements for the FMT reconstruction were approximately equal to the ratio between the light intensities measured at the emission and excitation wavelengths. Image reconstruction was carried out for a partial region of the phantom including the targets, with a height of 2 cm.

Figure 2 shows the fluorophore distribution reconstructed using the Tikhonov, re-L1-NCG, and ASDR-L1 methods, respectively. Figures 2(a) and 2(b) denote the cross section and 3D rendering of the true fluorophore distribution, where the two identical fluorophores with an EED of 0.6 cm were located at (0 cm, 1.0 cm, 2.4 cm) and (0 cm, 0 cm, 2.4 cm), respectively. In other words, one fluorophore was located very close to the edge of the phantom, while the other one was located at the center of the phantom. As can be seen from Figs. 2(d), 2(f) and 2(h), both the L1 regularization and the reweighted L1 regularization can reconstruct the fluorophore distribution with less noise compared with the Tikhonov regularization. Figures 2(c) and 2(d) show that the Tikhonov method results in an over-smoothing effect, i.e. the volume of the reconstructed targets is much larger than the true volume. And the reconstructed value of the center fluorophore is smaller than that of the edge fluorophore when Tikhonov is employed. As shown in Figs. 2(e) and 2(f), although re-L1-NCG as an L1 regularization method can shrink the reconstructed fluorophores, the edge and center fluorophores are reconstructed with different intensities. Moreover, neither Tikhonov nor re-L1-NCG can clearly separate the two identical fluorophores located at different depths. However, as shown in Figs. 2(g) and 2(h), ASDR-L1 as a reweighted L1-regularization method can reconstruct the fluorophore distribution with higher spatial resolution. In addition, the two identical fluorophores reconstructed using ASDR-L1 have similar intensities.

Fig. 2 Reconstruction results of the physical phantom experiments with an EED of 0.6 cm between two fluorophores. (a) Cross section of the true fluorophore distribution at the Z = 2.6 cm plane, (b) 3D rendering of true fluorescent targets, where the red circle denotes the position of the cross section. Slice (c) and stereo (d) are the reconstruction results obtained from Tikhonov. Slice (e) and stereo (f) correspond to re-L1-NCG. Slice (g) and stereo (h) correspond to ASDR-L1. Each image was normalized by its maximum of the reconstruction results.

Download Full Size | PDF

In order to verify the proposed algorithm when the fluorophores have a smaller EED, additional experiments were conducted. The two fluorophores had an EED of 0.4 cm and were located at (0 cm, 0.8 cm, 2.4 cm) and (0 cm, 0 cm, 2.4 cm), respectively. The results are shown in Fig. 3. As can be observed, re-L1-NCG and ASDR-L1 obtain the reconstruction results with less noise compared with Tikhonov. And by comparing Figs. 2 and 3, it can be seen that neither Tikhonov nor re-L1-NCG can obtain the separated fluorophores, when the distance between the edge and center fluorophores is reduced. However, ASDR-L1 can obtain the reconstruction results with higher spatial resolution and more consistent intensities at different depths, even if the EED of the two fluorophores decreases.

Fig. 3 Reconstruction results of the physical phantom experiment with an EED of 0.4 cm between two fluorophores. (a) Cross section of the true fluorophore distribution at the Z = 2.6 cm plane, (b) 3D rendering of true fluorescent targets, where the red circle denotes the position of the cross section. Slice (c) and stereo (d) are the reconstruction results obtained from Tikhonov. Slice (e) and stereo (f) correspond to re-L1-NCG. Slice (g) and stereo (h) correspond to ASDR-L1. Each image was normalized by its maximum of the reconstruction results.

Download Full Size | PDF

In order to better analyze the performance of the reconstruction methods in terms of spatial resolution and relative quantitation, the profiles across the two targets at X = 0 (i.e., along the Y axis) in the cross sections of Figs. 2 and 3 are shown in Figs. 4(a) and 4(b), respectively. As evident in Fig. 4, compared with the profiles corresponding to Tikhonov and re-L1-NCG, the FWHM of the profiles corresponding to ASDR-L1 is smaller and the valley value of the corresponding profile does not increase apparently when the EED decreases. This indicates that ASDR-L1 can obtain higher spatial resolution. For the profile corresponding to Tikhonov in Fig. 4(a), the peak value of the center target is about 30% of that of the edge target. This percentage is about 55% for re-L1-NCG in Fig. 4(a). This means that Tikhonov and re-L1-NCG have lower relative quantitation performance when two identical fluorophores are located at different depths. However, for the profiles corresponding to ASDR-L1 in Fig. 4, the reconstructed peak value of the center fluorophore is nearly equal to that of the edge fluorophore, which demonstrate the higher performance of ASDR-L1 in term of relative quantitation.

Fig. 4 Intensity profiles along the Y axis in the cross sections of (a) Figs. 2 and (b) 3.

Download Full Size | PDF

In order to further compare the relative quantitation performance of different reconstruction methods, the quantitative metrics calculated using Eq. (9) for the profiles in Fig. 4 are listed in Table 1. As shown, for the two cases, the RQ values of Tikhonov and re-L1-NCG are much larger than that of ASDR-L1, which indicates that ASDR-L1 has better performance in term of relative quantitation.

Table 1. Relative quantitation indexes calculated from Eq. (9) for the physical phantom experiments.

View Table | View all tables in this article

3.2 Mouse experiments

In order to further study the potential of the proposed method in practical applications, in vivo animal experiments were conducted under the protocol approved by the Institutional Animal Care and Use Committee of Tsinghua University. Two transparent glass tubes (0.5 cm in outer diameter, 1.0 cm in length) filled with 1.3 μM × 20 μL ICG were employed as the fluorescent targets. The two tubes were implanted in to the abdomen of an 8-week BALB/c nude mouse through surgical operation. Because the reconstruction region of interest is the abdomen, homogeneous model is adopted here. For the FMT reconstruction, the abdomen of nude mouse had an absorption coefficient of 2.38 × 10⁻² cm⁻¹ and reduced scattering coefficient of 10.3 cm⁻¹, calculated according to Ref [26]. The mouse experiments were conducted using the previously developed hybrid full-angle free-space FMT and X-ray micro-cone-beam computed tomography (CT) system [27], where the FMT module had the same setup as that used in the physical experiments. The CT images were obtained to examine the reconstruction accuracy of the proposed method. The mouse was anesthetized by isoflurane-oxygen gas mixture during the surgical operation and imaging process. 36 white-light images were collected to recover the 3D surface of the mouse abdomen [28]. The detector FOV for each projection was about 180° and 36 projections were collected during the 360 degree rotation.

The fusion result of the white-light and fluorescence images is shown in Fig. 5(a), where the two red dotted lines determinate the reconstruction range. The sensitivity matrix W obtained using the KA method had the size of 7097 × 12119, where the 3D geometry of the abdomen shown in Fig. 5(b) was discretized with a voxel size of 0.08 × 0.08 × 0.08 cm³ and the detector spatial sampling was set to 0.18 cm.

Fig. 5 Mouse experiments with two implanted fluorescence targets. (a) Fusion of the white-light and fluorescence images. Region between the two red dotted lines is used for reconstruction. (b) 3D solid geometry of the abdomen.

Download Full Size | PDF

Figure 6 shows the fluorophore distribution reconstructed using the Tikhonov, re-L1-NCG and ASDR-L1 methods, respectively. As shown, the fluorophores reconstructed using Tikhonov have too large volumes and much noise. Although re-L1-NCG as an L1 regularization method can obtain reconstruction result with higher spatial resolution than Tikhonov, the reconstructed intensity of the right fluorophore is lower than that of the left fluorophore, which indicates the poorer relative quantitation capability. Compared with the Tikhonov and re-L1-NCG methods, ASDR-L1 can obtain the reconstruction results with higher spatial resolution and less noise. Furthermore, the reconstructed intensities of the two identical fluorophores are similar, by means of ASDR-L1. To illustrate the reconstruction accuracy of the algorithms, the FMT slice images reconstructed using different methods are overlaid onto the corresponding CT slice. It can be observed that the proposed ASDR-L1 algorithm has fairly good reconstruction accuracy. It should be noted that the reconstructed location error is acceptable, because the errors may have been caused by the approximated forward model, the inaccuracy optical properties and other geometrical approximation, such as the rough point-to-point mapping between the captured two dimensional (2D) CCD images and the light flux on the 3D mouse surface.

Fig. 6 Reconstruction results of the mouse experiments. Columns (a)-(c) represent the slice and stereo results reconstructed using Tikhonov, re-L1-NCG and ASDR-L1, respectively. Each slice is composed of the reconstructed FMT slice fused with CT slice. Each image was normalized by its maximum of the reconstruction results.

Download Full Size | PDF

In order to better compare the performance of different reconstruction methods, the profiles along the white dotted lines in Fig. 6 are shown in Fig. 7. It can be seen that the reconstruction results obtained using ASDR-L1 and re-L1-NCG have smaller FWHM compared with Tikhonov. For the Tikhonov and re-L1-NCG methods, the reconstructed intensities of the two identical fluorophores are significantly different. For the profile corresponding to ASDR-L1, the peak value of the right fluorophore is nearly 95% of that of the left fluorophore. This percentage is about 58% and 45% for Tikhonov and re-L1-NCG, respectively. Thus, ASDR-L1 has better performance in terms relative quantitation and spatial resolution, which is consistent with the results of physical phantom experiments.

Fig. 7 Intensity profiles along the white dotted lines in the cross sections of Fig. 6.

Download Full Size | PDF

The relative quantitation indexes for the profiles in Fig. 7 are listed in Table 2, so as to quantitatively demonstrate the performance of the proposed ASDR-L1 algorithm. As shown in Table 2, ASDR-L1 has much smaller RQ value than Tikhonov and re-L1-NCG, which demonstrates the improved performance of ASDR-L1 in term of relative quantitation.

Table 2. Relative quantitation indexes calculated from Eq. (9) for the mouse experiments.

View Table | View all tables in this article

4. Discussion and conclusion

The ill-posedness of the FMT inverse problem causes relatively low spatial resolution in the reconstruction results, especially when L2 regularization methods are used. In addition, the nonlinear depth sensitivity of the measurements causes the inconsistent reconstructed intensities of the fluorophores with identical concentration. In this work, ASDR-L1 as an adaptively reweighted L1 regularization algorithm is proposed to compensate for the spatial dependence of the intensity and resolution. The ASDR-L1 algorithm is composed of the adaptive support estimate and the iteratively updated weights. The adaptive support estimate is updated based on the threshold strategy in the reweighted RLI loop of ASDR-L1. And the support estimate can be regarded as the permission region for the reconstructed fluorescent targets. In the corresponding support estimate, the weights for the reweighted L1 regularization sub-problem are iteratively updated as a function of the preceding iterative solution. Actually, for FMT reconstruction, the iteratively updated weights in ASDR-L1 mean different regularization parameters at different spatial locations. Thus to some extent, ASDR-L1 can be considered as one kind of spatially variant regularization methods [19]. However, unlike the method proposed in Ref [19], the spatially variant regularization parameters in ASDR-L1 are updated adaptively and independent of radius information. Furthermore, it should be noted that ASDR-L1 also can be seen one kind of non-convex regularization algorithms which are recently studied in FMT to enhance the sparsity [12,13].

In order to compare the performance among L2, L1 and reweighted L1 regularizations, the Tikhonov and re-L1-NCG methods are employed in this work. Physical phantom experiments with two fluorophores located at different depths were conducted to compare the performances. The results demonstrate that ASDR-L1 can obtain reconstruction results with consistent intensities at different depths and high spatial resolution. In addition, the in vivo mouse experiments were carried out to further demonstrate the performance of ASDR-L1. Unlike the physical phantom experiments, in the mouse experiments, the depths of the two targets have no obvious differences. Although ASDR-L1 has better relative quantitation performance than Tikhonov and re-L1-NCG, the reconstructed values of the two identical fluorophores are still slightly different. This may be caused by the irregular 3D surface of the mouse abdomen, the inaccuracy of the optical properties, the approximated forward model and other factors which will be studied in the future.

In order to further improve the localization accuracy and relative quantitation for FMT reconstruction, more accurate models to describe photon propagation in biological tissues and free space are necessary [20,29]. Except for DE, the forward model based on RTE or simplified spherical harmonics approximation can be adopted, while the heterogeneous model can be utilized for more accurate solution. An effective method can be employed to accurately reconstruct photon flux on a body surface of arbitrary shape from 2D photographic images measured by a CCD camera [29]. The ASDR-L1 algorithm can potentially be utilized in FMT reconstruction with these improved models.

In conclusion, ASDR-L1 as one kind of reweighted L1-regularization methods has been proposed to solve the non-convex regularization problem for FMT. And the reconstruction results of physical phantom and mouse experiments demonstrate the improved performance of ASDR-L1 in terms of relative quantitation and spatial resolution. Future works will be focused on the in vivo small animal experiments with probe-marked tumor models. And in order to realize an effective iterative reconstruction algorithm, a robust termination criterion will be studied in our future work.

Acknowledgments

This work is supported by the National Basic Research Program of China (973) under Grant No. 2011CB707701; the National Natural Science Foundation of China under Grant No. 81227901, 81271617, 61322101, 61361160418, 61401246; the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114; the China Postdoctoral Science Foundation under Grant No. 2014M550073; Fei Liu is supported in part by the Postdoctoral Fellowship of Tsinghua-Peking Center for Life Sciences.

References and links

1. C. Darne, Y. Lu, and E. M. Sevick-Muraca, “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. 59(1), R1–R64 (2014). [CrossRef] [PubMed]

2. S. R. Cherry, “In vivo molecular and genomic imaging: new challenges for imaging physics,” Phys. Med. Biol. 49(3), R13–R48 (2004). [CrossRef] [PubMed]

3. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (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. J. Haller, D. Hyde, N. Deliolanis, R. de Kleine, M. Niedre, and V. Ntziachristos, “Visualization of pulmonary inflammation using noninvasive fluorescence molecular imaging,” J. Appl. Physiol. 104(3), 795–802 (2008). [CrossRef] [PubMed]

7. J. Shi, B. Zhang, F. Liu, J. Luo, and J. Bai, “Efficient L1 regularization-based reconstruction for fluorescent molecular tomography using restarted nonlinear conjugate gradient,” Opt. Lett. 38(18), 3696–3699 (2013). [CrossRef] [PubMed]

8. J. Shi, F. Liu, G. Zhang, J. Luo, and J. Bai, “Enhanced spatial resolution in fluorescence molecular tomography using restarted L1-regularized nonlinear conjugate gradient algorithm,” J. Biomed. Opt. 19(4), 046018 (2014). [CrossRef] [PubMed]

9. P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46(10), 1679–1685 (2007). [CrossRef] [PubMed]

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

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

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

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

14. E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted L1 minimization,” Journal of Fourier Analysis and Applications 14(5–6), 877–905 (2008). [CrossRef]

15. H. Mansour and O. Yilmaz, “Support driven reweighted L1 minimization,” IEEE International Conference on Speech and Signal Processing (ICASSP), IEEE Press: Piscataway NJ, 2012: 3309–3312.

16. Y. Wang and W. Yin, “Sparse signal reconstruction via iterative support detection,” SIAM Journal on Imaging Sciences 3(3), 462–491 (2010). [CrossRef]

17. D. Wipf and S. Nagarajan, “Iterative reweighted and methods for finding sparse solutions,” IEEE J. Sel. Top. Sig. Processing 4(2), 317–329 (2010). [CrossRef]

18. C. B. Shaw and P. K. Yalavarthy, “Performance evaluation of typical approximation algorithms for nonconvex Lp-minimization in diffuse optical tomography,” J. Opt. Soc. Am. A 31(4), 852–862 (2014). [CrossRef]

19. B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Osterberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38(13), 2950–2961 (1999). [CrossRef] [PubMed]

20. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), 41–93 (1999). [CrossRef]

21. A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging 24(10), 1377–1386 (2005). [CrossRef] [PubMed]

22. J. Ripoll, M. Nieto-Vesperinas, R. Weissleder, and V. Ntziachristos, “Fast analytical approximation for arbitrary geometries in diffuse optical tomography,” Opt. Lett. 27(7), 527–529 (2002). [CrossRef] [PubMed]

23. J. Shi, F. Liu, J. Luo, and J. Bai, “Depth compensation in fluorescence molecular tomography using an adaptive support driven reweighted L1-minimization algorithm,” Proc. SPIE 9216, 921603 (2014). [CrossRef]

24. F. Liu, M. Li, B. Zhang, J. Luo, and J. Bai, “Weighted depth compensation algorithm for fluorescence molecular tomography reconstruction,” Appl. Opt. 51(36), 8883–8892 (2012). [CrossRef] [PubMed]

25. F. Liu, X. Liu, D. Wang, B. Zhang, and J. Bai, “A parallel excitation based fluorescence molecular tomography system for whole-body simultaneous imaging of small animals,” Ann. Biomed. Eng. 38(11), 3440–3448 (2010). [CrossRef] [PubMed]

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

27. X. Guo, X. Liu, X. Wang, F. Tian, F. Liu, B. Zhang, G. Hu, and J. Bai, “A combined fluorescence and microcomputed tomography system for small animal imaging,” IEEE Trans. Biomed. Eng. 57(12), 2876–2883 (2010). [CrossRef] [PubMed]

28. D. Wang, X. Liu, Y. Chen, and J. Bai, “In-vivo fluorescence molecular tomography based on optimal small animal surface reconstruction,” Chin. Opt. Lett. 8(1), 82–85 (2010). [CrossRef]

29. X. Chen, X. Gao, D. Chen, X. Ma, X. Zhao, M. Shen, X. Li, X. Qu, J. Liang, J. Ripoll, and J. Tian, “3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images,” Opt. Express 18(19), 19876–19893 (2010). [CrossRef] [PubMed]

Reconstruction Methods	RQ for Fig. 4(a)	RQ for Fig. 4(b)
Tikhonov	0.2622	0.1734
Re-L1-NCG	0.1472	0.1203
ASDR-L1	0.0020	0.0035

Reconstruction Methods	RQ for Fig. 7
Tikhonov	0.1309
Re-L1-NCG	0.1935
ASDR-L1	0.0139

Reconstruction Methods	RQ for Fig. 4(a)	RQ for Fig. 4(b)
Tikhonov	0.2622	0.1734
Re-L1-NCG	0.1472	0.1203
ASDR-L1	0.0020	0.0035

Reconstruction Methods	RQ for Fig. 7
Tikhonov	0.1309
Re-L1-NCG	0.1935
ASDR-L1	0.0139

An adaptive support driven reweighted L1-regularization algorithm for fluorescence molecular tomography

Abstract

1. Introduction

2. Methods

2.1 Forward model for FMT

2.2 FMT reconstruction based on reweighted L1 regularization

2.3 Quantitative metrics for FMT reconstruction

3. Experiments and results

3.1 Physical phantom experiments

3.2 Mouse experiments

4. Discussion and conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Tables (2)

Equations (9)

Biomedical Optics Express