Fluorescence molecular tomography reconstruction via discrete cosine transform-based regularization

Junwei Shi; Fei Liu; Jiulou Zhang; Jianwen Luo; Jing Bai

doi:10.1117/1.JBO.20.5.055004

13 May 2015 Fluorescence molecular tomography reconstruction via discrete cosine transform-based regularization

Junwei Shi, Fei Liu, Jiulou Zhang, Jianwen Luo, Jing Bai

Author Affiliations +

Journal of Biomedical Optics, Vol. 20, Issue 5, 055004 (May 2015). https://doi.org/10.1117/1.JBO.20.5.055004

Abstract

Fluorescence molecular tomography (FMT) as a noninvasive imaging modality has been widely used for biomedical preclinical applications. However, FMT reconstruction suffers from severe ill-posedness, especially when a limited number of projections are used. In order to improve the quality of FMT reconstruction results, a discrete cosine transform (DCT) based reweighted L1-norm regularization algorithm is proposed. In each iteration of the reconstruction process, different reweighted regularization parameters are adaptively assigned according to the values of DCT coefficients to suppress the reconstruction noise. In addition, the permission region of the reconstructed fluorophores is adaptively constructed to increase the convergence speed. In order to evaluate the performance of the proposed algorithm, physical phantom and in vivo mouse experiments with a limited number of projections are carried out. For comparison, different L1-norm regularization strategies are employed. By quantifying the signal-to-noise ratio (SNR) of the reconstruction results in the phantom and in vivo mouse experiments with four projections, the proposed DCT-based reweighted L1-norm regularization shows higher SNR than other L1-norm regularizations employed in this work.

1. Introduction

Fluorescence molecular tomography (FMT) is a noninvasive biomedical imaging modality that exploits the specificity of fluorescent biomarkers to monitor three-dimensional (3-D) cell and molecular activities inside biological tissues.¹ As a rapidly growing technique, FMT has already been applied to clinical and preclinical studies.²^,³ However, owing to the high scattering of light propagating through tissues, the FMT inverse problem is inherently ill-posed and quite susceptible to measurement data noise and model errors. This makes accurate and robust FMT reconstruction a challenge, especially for FMT reconstruction with only a limited number of projections. In order to obtain high-quality FMT reconstruction results, additional prior information is usually incorporated in the FMT inverse problem through some regularization techniques, including L2-norm,⁴ L1-norm,⁵ Lp norm ( $0 < p < 1$ )⁶, total variation (TV) regularizations, etc.⁷ Although L2-norm regularization is simple to implement and can be efficiently solved, it results in blurred images with low spatial resolution and produces noisy background when the inverse problem is severely ill posed.⁴ Based on the fact that fluorophore distribution is usually sparse (such as the early-stage tumors tagged with fluorescent probes), L1-norm regularization as a convex relaxation of L0-norm regularization has been extensively studied and found to be able to yield results with high-spatial resolution.⁵ In addition, as the shape-based reconstruction method, the sparsity-promoting TV regularization⁷ and spherical harmonics parameterization⁸ can be employed to preserve edges in images for piecewise-constant reconstruction. Recently, as a nonconvex approximation to L0-norm regularization, the reweighted L1-norm scheme and Lp-norm regularization have been reported to outperform the regular L1-norm regularization.⁶^,⁹ However, some consistent efforts still need to be made to improve the reconstruction quality for FMT, especially when only a limited number of projections are used.

In this paper, a discrete cosine transform (DCT) based reweighted L1-norm regularization (DCT-re-L1) algorithm is proposed to improve the quality of reconstruction results for the ill-posed FMT inverse problem. As an orthogonal transform, DCT is an effective tool for image compression because of its energy compaction property (i.e., most of the signal energy focuses on a few large DCT coefficients).¹⁰ In view of the high degree of space correlation, the spatial distribution of fluorescent targets can be effectively compressed by 3-D DCT. Thus, the sparsity character in the 3-D DCT domain can be seen as a kind of prior information to guarantee the uniqueness and stability of the FMT inverse problem. In order to suppress the reconstruction noise, the high-frequency coefficients in the DCT domain should be penalized more heavily because they usually correspond to noise. The low-frequency coefficients in the DCT domain should be assigned with small penalties, because most of the image information is retained in the low-frequency coefficients.⁹ In order to suppress the reconstruction noise in each iteration, an iteratively updated scheme is proposed to adaptively assign different L1-regularization parameters to different DCT coefficients. In addition, we propose a scheme to update the permission region (without additional anatomical information) based on the non-negative region of the results in the previous iterations, which can speed up the convergence and computational process and reduce the memory consumption.

This paper is organized as follows. In Sec. 2, the mathematical framework of FMT and the proposed algorithm are presented. In Sec. 3, the physical phantom and in vivo mouse experiments with different numbers of FMT projections are conducted to evaluate the performance of the proposed algorithm. The results are discussed and this paper is concluded in Sec. 4.

2. Methods

2.1.

Forward Model for FMT

The propagation of near-infrared light in highly scattering tissues is modeled by the coupled diffusion equations (DE) with Robin-type boundary condition, when the sources and detectors are separated by more than a few mean-free-path lengths.¹¹ In this paper, the normalized Born approximation to DE is adopted to linearize the nonlinear FMT problem. This normalization scheme can make the FMT forward model insensitive to theoretical errors about boundary conditions and tissue optical properties. Through Born approximation to DE, 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 calculated as follows:¹²

Eq. (1)

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

where

ϕ_{em} (r_{d}, r_{s})

and

ϕ_{ex} (r_{d}, r_{s})

denote the average intensity at the emission wavelength

λ_{2}

and excitation wavelength

λ_{1}

, respectively.

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

is the theoretically calculated photon field at the excitation wavelength

λ_{1}

induced at position

r

by a source at position

r_{s}

.

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

is the Green’s function describing photon propagation from point

r

to the detector point

r_{d}

at the emission wavelength

λ_{2}

.

D^{f}

is the diffusion coefficient of the tissue at wavelength

λ_{2}

.

υ

is the speed of light in the tissue.

S_{0}

is a calibration factor accounting for the system gain and attenuation. After the image domain is discretized, the FMT forward model can be analytically solved by the Kirchhoff approximation method,¹³ and the FMT problem can be formulated as the following linear matrix equation:

Eq. (2)

Φ_{f} = W X,

where

W \in R^{m \times n}

is the sensitive matrix mapping the vectorized fluorophore concentration

X \in R^{n \times 1}

into the fluorescence measurement vector

Φ_{f} \in R^{m \times 1}

.

2.2.

DCT-Regularized Reconstruction Algorithm

In order to obtain high-quality reconstruction results from the ill-posed Eq. (2), sparsity-promoting regularization based on DCT is employed and the corresponding FMT inverse problem can be formulated as the following optimization function:

Eq. (3)

\underset{X \geq 0}{argmin} ({‖ W X - Φ_{f} ‖}_{2}^{2} + Λ^{k} {‖ D X ‖}_{1}) .

Equation (3) contains two parts: the data fidelity term which is used to match the solution with the measurements and the regularization term which represents a kind of prior penalty. Matrix $D \in R^{L \times n}$ denotes the real 3-D DCT matrix for the fluorophore-concentration vector $X$ , where $L$ is the number of frequency coefficients in the 3-D DCT domain. $Λ^{k}$ is the diagonal regularization-parameter matrix updated in the $k$ ’th iteration, with $λ_{i}^{k} (i = 1, 2, \dots, L)$ on the diagonal and zero elsewhere. In order to roughly level off the differences in detection sensitivities, Eq. (3) can be rewritten as follows:⁵

Eq. (4)

\underset{X \geq 0}{argmin} ({‖ W_{new} X - Φ_{fnew} ‖}_{2}^{2} + Λ^{k} {‖ D X ‖}_{1}),

Eq. (5)

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

Eq. (6)

W_{new} = W W_{nor}; Φ_{fnew} = Φ_{f} / \max (Φ_{f}),

where

W_{i}

denotes the

i

’th column of the sensitive matrix

W

. In order to solve the DCT regularized Eq. (4), the proposed DCT-re-L1 algorithm based on an adaptively updated regularization parameter and permission region is proposed. As is similar to our previously proposed restarted L1 regularization-based nonlinear conjugate gradient (re-L1-NG) algorithm,⁵ in each iteration of DCT-re-L1, the optimization function is solved by the L1-regularized nonlinear conjugate gradient (L1-NCG) method with a backtracking line search.⁵ The main steps of the proposed DCT-re-L1 algorithm are summarized in Algorithm 1:

Algorithm 1

DCT-based reweighted L1-norm regularization.

Step 1: Initialization: Set

k_{inner} = 100

as the maximum number of inner iterations for L1-NCG,

σ = 0.0001

as the threshold for the termination condition of DCT-re-L1.

k = 1

. In L1-NCG, other parameters such as the initial value, original searching direction, and two coefficients about the backtracking line search can be found in Ref. 5.

Step 2: For Eq. (4) without a non-negative constraint, compute the solution

X_{L 1 - NCG}^{1}

using L1-NCG with

k_{inner}

iterations, where the regularization-parameter matrix

Λ^{1}

in Eq. (4) is initialized as an identity matrix. Then set the negative value in

X_{L 1 - NCG}^{1}

to be zero.

Step 3: The permission region is initialized as

{per}^{1} = {1,2, \dots, n}

. The first nonpositive region is set as

{neg}^{1} = {indices of the zero value in X_{L 1 - NCG}^{1}}

. The normalized norm of the residual is calculated as

{RMS}^{1} = norm (W_{new} X_{L 1 - NCG}^{1} - Φ_{fnew}) / norm (Φ_{fnew})

.

Step 4: In the matrices

W_{new}

and

D

, only keep the columns whose indices are incorporated in

{per}^{k}

and obtain the new matrices

W_{new}^{per}

and

D^{per}

. Then the new optimization function based on the permission region turns into:

Eq. (7)

\underset{X \geq 0}{argmin} ({‖ W_{new}^{per} X^{per} - Φ_{fnew} ‖}_{2}^{2} + Λ^{k + 1} {‖ D^{per} X^{per} ‖}_{1})

where

X^{per} \in R^{p \times 1}

denotes the reconstruction-result vector in the permission region and

p

is the number of indices in

{per}^{k}

.

Step 5: The initial value for Eq. (7) is set as

X_{0}^{per} = X_{L 1 - NCG}^{k} ({per}^{k})

. Calculate the DCT-coefficient vector

E = D^{per} X_{0}^{per}

(

E \in R^{P \times 1}

).

Max = \max (| E |)

. Generate the updated regularization-parameter matrix

Λ^{k + 1}

, with

λ_{i}^{k + 1} = g (e_{i}, Max)

on the diagonal and zero elsewhere.

Eq. (8)

g (e_{i}, Max) = {\begin{matrix} Max / (0.01 \times Max) & if | e_{i} | \leq 0.01 \times Max \\ Max / | e_{i} | & if | e_{i} | > 0.01 \times Max \end{matrix}

Step 6: Solve Eq. (7) for

X^{per} \in R^{p \times 1}

using L1-NCG with

k_{inner}

iterations. Then obtain

X_{L 1 - NCG}^{k + 1} \in R^{n \times 1}

where the values corresponding to the indices in

{per}^{k}

are

X^{per} \in R^{p \times 1}

, and the other values are zero.

Step 7: Set the negative value in

X_{L 1 - NCG}^{k + 1}

to be zero. Construct the nonpositive region as:

{neg}^{k + 1} = {indices of the zero values in X_{L 1 - NCG}^{k + 1}}

. Update the normalized residual norm as:

Eq. (9)

{RMS}^{k + 1} = norm (W_{new} X_{L 1 - NCG}^{k + 1} - Φ_{fnew}) / norm (Φ_{fnew}) .

Construct the corresponding permission region as

Eq. (10)

{per}^{k + 1} = \bar{{neg}^{k} \cap {neg}^{k + 1}}

Step 8:

k = k + 1

. If

| {RMS}^{k} - {RMS}^{k + 1} | / {RMS}^{k} < σ

or

k > k_{\max}

(maximum number of outer iterations), stop the iteration. Otherwise, go to Step 4.

The first three steps are used to initialize the permission region and the residual by setting an identity matrix as the initialized regularization-parameter matrix $Λ^{1}$ . In the first iteration, all the DCT coefficients are penalized equally because there is no prior information about the fluorophore distribution. Steps 4 to 8 constitute the main iterative loop, where the regularization-parameter matrix $Λ^{k}$ and the permission region ${per}^{k}$ are adaptively updated. In order to suppress the reconstruction noise, the high-frequency information corresponding to the small DCT coefficients are penalized heavily, which is the basis by which to construct $λ_{i}^{k}$ as the reciprocal value of the previous DCT coefficient. In order to avoid generating infinite $λ_{i}^{k}$ when the corresponding DCT coefficient $e_{i}$ is near zero, parameter $λ_{i}^{k}$ is set to be 100 [i.e., $Max / (0.01 \times Max)$ ] when $| e_{i} | \leq 0.01 \times Max$ . As shown in Step 5, Max denotes the maximum value of the calculated DCT coefficients in each iteration. Owing to the non-negativity constraint of Eq. (4), the negative values in $X_{L 1 - NCG}^{k}$ are set to be zero, and this zero region ${neg}^{k}$ can provide useful information about the permission region. As shown in Eq. (10), the permission region ${per}^{k + 1}$ is iteratively updated as the complement of set ${neg}^{k} \cap {neg}^{k + 1}$ , which means that if a reconstruction value is negative in two consecutive iterations, its index will not belong to the permission region. Because the size of the permission region is smaller than that of the entire reconstruction region, the dimension of the sensitive matrix in Eq. (7) scales down compared with Eq. (4), which reduces the memory consumption and increases the computational speed. The algorithm is terminated when the relative change of the residual norm between two consecutive iterations is less than $σ$ .

In order to demonstrate the performance differences between the regular L1-norm regularization and the DCT-based reweighted L1-norm regularization, Eq. (4) is compared with the following optimization function:⁵

Eq. (11)

\underset{X \geq 0}{argmin} ({‖ W_{new} X - Φ_{fnew} ‖}_{2}^{2} + λ {‖ X ‖}_{1}) .

In order to demonstrate the advantage of the proposed DCT-based reweighted strategy, Eq. (4) is compared with the following DCT regularized optimization function:

Eq. (12)

\underset{X \geq 0}{argmin} ({‖ W_{new} X - Φ_{fnew} ‖}_{2}^{2} + λ {‖ D X ‖}_{1}) .

Unlike Eq. (4), a constant regularization parameter is assigned to all DCT coefficients in Eq. (12). In addition, in order to compare the performance between different reweighted L1-norm strategies, Eq. (4) is compared with the classic reweighted L1-norm regularization function as follows:⁹^,¹⁴

Eq. (13)

\underset{X \geq 0}{argmin} ({‖ W_{new} X - Φ_{fnew} ‖}_{2}^{2} + λ {‖ R^{k} X ‖}_{1}),

where

R^{k}

is the diagonal weight matrix in the

k

’th iteration, with

r_{i}^{k} (i = 1, 2, \dots, L)

on the diagonal and zeros elsewhere. The weighed matrix

R^{k}

is iteratively updated as follows:

Eq. (14)

r_{i}^{k} = \frac{1}{| x_{i}^{k} | + α},

where

α

is a stability parameter that used to avoid the division by a zero-value component

x_{i}^{k}

. In this paper,

α

is empirically set to be 0.01.⁹ In order to better focus on the comparison between different L1-norm regularization strategies, the solving algorithms for Eqs. (11) to (13) are the same as Algorithm 1 except for the difference in Step 5. When solving Eqs. (11) to (13), the corresponding regularization parameters are manually optimized to be 30. All the reconstruction algorithms are implemented in MATLAB 2008a (MathWorks, Natick, Massachusetts) on a personal computer with a i7-2600 CPU (3.40 GHz, Intel, Santa Clara, California) and 8 GB RAM.

3. Experiments and Results

3.1.

Physical Phantom Experiments

In order to evaluate the performance of the proposed algorithm, physical phantom experiments were performed with the previously developed dual-modality FMT/X-ray computed tomography (CT) imaging system.¹⁵ The schematic diagram of the system is shown in Fig. 1, where the custom-made free-space FMT and micro-cone-beam-CT subsystems were orthogonally placed on an optical bench.

Fig. 1

Schematic top view of the dual-modality fluorescence molecular tomography (FMT)/computed tomography (CT) system.

In the physical phantom experiments, a transparent glass cylinder with a diameter of 3 cm and height of 6 cm filled with 1% intralipid was employed as the phantom. The absorption and reduced scattering coefficients were 0.02 and $10 {cm}^{- 1}$ , respectively. The fluorescence targets in the phantom were two identical cylinders with diameters of 0.4 cm, each of which was filled with $10 μ L$ indocyanine green (ICG) with a concentration of $1.3 μ M$ . They were immersed in the phantom with an edge-to-edge distance of 0.8 cm, and their centers were positioned at (0.22 cm, $- 0.55 cm$ , 2.65 cm) and ( $- 0.21 cm$ , 0.58 cm, 2.65 cm), respectively. Through the customized optical fiber, a point incident light was generated. In order to excite ICG, a bandpass filter with a center wavelength of 770 nm and full-width at half maximum (FWHM) of 10 nm was placed in front of a Xenon lamp (Asahi Spectra, Torrance, California). In order to collect the fluorescence-light measurement, a bandpass filter with a center wavelength of 840 nm and FWHM of 12 nm was placed in front of the highly-sensitive cooled charge coupled device (CCD) camera (iXon DU-897, Andor Technologies, Belfast, Northern Ireland). When collecting the excitation-light measurement, a neutral density filter (with an optical density of 2) in place of the fluorescence filter was placed in front of the CCD camera. The detector field of view corresponding to a point excitation source was about 180 deg and the detector spatial sampling was set to be 0.2 cm. The normalized Born measurements in Eq. (1) were calculated from the ratio between the measurements corresponding to the fluorescence and excitation light.

First, the performances of different L1-norm regularizations were compared based on the reconstruction results with 18 projections. The dimensions of the corresponding sensitive matrix $W$ were $2984 \times 6196$ , where the discretized mesh resolution of the phantom was set as $0.13 \times 0.13 \times 0.13 {cm}^{3}$ . The true geometry configuration of the phantom is shown in Figs. 2(a1) and 2(a2). The red circle in Fig. 2(a2) denotes the position of the reconstructed FMT/CT slice. The second column of Fig. 2 denotes the results of Eq. (4) (i.e., DCT-based reweighted L1-norm regularization function). The third column is obtained by solving Eq. (12) (i.e., DCT-based L1-norm regularization function without reweighted strategy). The fourth and fifth columns correspond to Eq. (11) (i.e., the regular L1-norm regularization function) and Eq. (13) (i.e., the classic reweighted L1-norm regularization function), respectively.

Fig. 2

Reconstruction results of the phantom experiments with 18 projections. (a1) Slice and (a2) three-dimensional (3-D) rendering of the true fluorescent targets. (b1) and (b2) Slice and 3-D results of DCT-based reweighted L1-norm regularization. (c1) and (c2) Slice and 3-D results of DCT-based L1-norm regularization without reweighted strategy. (d1) and (d2) Slice and 3-D results of regular L1-norm regularization. (e1) and (e2) Slice and 3-D results of classic L1-norm regularization. The red circle in (a2) denotes the position of the FMT/CT slice. (f) Each slice is composed of the reconstructed FMT result fused with corresponding CT slice. Each image is normalized by its maximum.

By overlaying the reconstructed FMT results onto the corresponding CT slices in Figs. 2(b1), 2(c1), 2(d1), and 2(e1), it can be seen that all the L1-norm regularizations can accurately locate the fluorescence targets. In order to further demonstrate the performances of different regularization techniques, the profiles along the white dotted lines in Fig. 2 are shown in Fig. 3. In Fig. 2, it can be seen that the lines in FMT cross-sections pass through the centers of the reconstructed and true targets at the same time. In Fig. 3, it can be seen that the profiles corresponding to all L1-norm regularizations have similar FWHMs. This means that these employed L1-norm regularizations can reconstruct the targets with similar spatial resolution. For all the profiles, the peak positions (corresponding to the center of the reconstructed target) match the true locations, which mean that all the employed L1-norm regularizations can obtain a high localization accuracy.

Fig. 3

Intensity profiles along the white dotted lines in the cross-sections of Fig. 2.

Although the cross-section results and profiles demonstrate similar performances of different methods, the comparison among the reconstructed 3-D images [Figs. 2(b2), 2(c2), 2(d2), and 2(e2)] demonstrates the performance differences of different L1-norm regularization strategies. The comparison between Figs. 2(b2) and 2(d2) validates the superior performance of the DCT-based reweighted L1-norm regularization over the regular L1-norm regularization in terms of noise suppressing. The comparison between Figs. 2(b2) and 2(c2) clarifies the fact that the better performance of the proposed DCT-re-L1 algorithm mainly relies on the proposed DCT-based reweighted strategy rather than the simple DCT-based regularization. Through the comparison between Figs. 2(d2) and 2(e2), it can be seen that the classic reweighted L1-norm regularization has superiority over the regular L1-norm regularization in terms of noise suppressing.

In order to further demonstrate the performance of the proposed DCT-re-L1 algorithm, limited-projection FMT reconstructions with nine, four, and three projections were carried out, respectively. Since the results obtained with Eq. (13) outperforms the results obtained with Eqs. (11) and (12) in terms of noise suppressing [Figs. 2(e2), 2(d2), and 2(c2)], only the DCT-based and classic reweighted L1-norm regularizations are compared. Figure 4 shows the results for different reweighted L1-norm regularizations with nine, four, and three projections, where the dimensions of the corresponding sensitive matrix $W$ were $1340 \times 6196$ , $678 \times 6196$ , and $495 \times 6196$ , respectively. The reconstructed slices shown in Figs. 4(b1) and 4(d1) demonstrate the good localization accuracy of the classic reweighted L1-norm regularization for limited-projection FMT with nine and four projections. However, Fig. 4(d2) shows that the 3-D image reconstructed by the classic reweighted L1-norm regularization is contaminated with noise when the number of FMT projections decreases to 4. That is partly because the corresponding sensitive matrix $W$ becomes severely underdetermined for the four-projection FMT reconstruction. In contrast, the results shown in Figs. 4(a1), 4(a2), 4(c1), and 4(c2) demonstrate that the DCT-based reweighted L1-norm regularization still shows good performances in terms of localization accuracy and noise suppressing, even for the four-projection FMT reconstruction. Compared with the 18-projection FMT, four-projection FMT only needs about one-fifth of the data-acquisition time to obtain high-quality reconstruction results. Additionally, because the number of measurements is reduced, the memory consumption is also reduced for the four-projection FMT. The reconstruction time of the DCT-based reweighted L1-norm regularization with 18 and 4 projections is 314 and 194 s, respectively. By contrast, the reconstruction time of the classic reweighted L1-norm regularization with 18 and 4 projections is 305 and 185 s, respectively. The 3-D DCT transform in Eq. (4) consumes more time in the reconstruction using the DCT-based regularization. Figures 4(e1), 4(e2), 4(f1), and 4(f2) show that the results of both the reweighted L1-norm regularizations have low localization accuracy and much noise when the number of FMT projections decreases to 3, which is caused by the insufficient fluorescence measurements used in the reconstruction process.

Fig. 4

Reconstruction results of the phantom experiments with nine, four, and three projections. Slices (a1, c1, e1) and 3-D (a2, c2, e2) denote the results reconstructed by the DCT-based reweighted L1-norm regularization. Slice (b1, d1, f1) and 3-D (b2, d2, f2) are the results reconstructed by the regular reweighted L1-norm regularization. Each slice is composed of the FMT result fused with CT slice. Each FMT image is normalized by its maximum.

In order to better analyze the performance of the DCT-based and classic reweighted L1-norm regularizations when the number of projections is reduced (i.e., from 18 to 3), the corresponding profiles along the white dotted lines in Figs. 2 and 4 are shown in Fig. 5. It can be seen that the width of the valleys decreases when the number of projections is reduced. This means that the smaller number of projections deteriorates the spatial resolution, especially when the number of projections is reduced to three [Fig. 5(d)]. As shown in Figs. 4(c1), 4(d1), 4(e1), and 4(f1), although the lines in FMT cross-sections pass through the centers of true targets, the centers of the reconstructed targets deviate these lines. This indicates that the centers of the reconstructed targets deviate from the true target centers.

Fig. 5

Intensity profiles along the white dotted lines in the cross-sections reconstructed by the DCT-based and classic reweighted L1-norm regularizations with different number of projections. The profiles in (a) correspond to Figs. 2(b1) and 2(e1) (18 projections). The profiles in (b) correspond to Figs. 4(a1) and 4(b1) (nine projections). The profiles in (c) correspond to Figs. 4(c1) and 4(d1) (four projections). The profiles in (d) correspond to Figs. 4(e1) and 4(f1) (three projections).

In order to quantitatively compare the localization accuracy, the center-localization errors corresponding to different reweighted L1-norm regularizations with different numbers of projections are listed in Table 1. As shown, the localization accuracy declines when the number of projections is reduced. When the number of projections decreases from nine to four, the localization errors have a dramatic increase.

Table 1

Center-localization errors (cm) of the results reconstructed by different reweighted L1-norm regularizations with different numbers of projections for the phantom experiments.

Number of projections	Classic reweighted L1 regularization	DCT-based reweighted L1 regularization
18	0.026	0.023
9	0.031	0.026
4	0.092	0.061
3	0.127	0.123

In order to quantitatively analyze the performance of different reweighted L1-norm regularizations in terms of noise suppressing, the signal-to-noise ratio (SNR) was calculated as

Eq. (15)

SNR = 10 \lg \frac{{‖ X (r_{target}) ‖}_{2}}{{‖ X (r_{back}) ‖}_{2}},

where

X (r_{target})

and

X (r_{back})

denote the reconstructed fluorophore distribution at target and background locations, respectively. The target and background location was defined by the true fluorophore distribution shown in Figs. 2(a1) and 2(a2). The corresponding SNR values for FMT reconstruction with 18, 9, 4, and 3 projections are calculated and listed in Table 2. As shown, the results reconstructed by the DCT-based reweighted L1-norm regularization have higher SNRs than those reconstructed by the classic reweighted L1-norm regularization, especially when four projections are used for FMT reconstruction. In addition, Table 2 shows that there is a big drop of the SNR values corresponding to both DCT-based and classic reweighted L1-norm regularizations when the number of projections decreases from nine to four. In addition to the reconstructed noise, such a drop is partly caused by the slight localization errors (shown in Table 1) and the reduced spatial resolution (shown in Fig. 5), when the number of projections used in FMT reconstruction is reduced to four.

Table 2

Signal-to-noise ratio (SNR) of the results reconstructed by different reweighted L1-norm regularizations with different numbers of projections for the phantom experiments.

Number of projections	Classic reweighted L1 regularization	DCT-based reweighted L1 regularization
18	11.5	12.1
9	9.7	11.2
4	3.5	8.4
3	$- 1.6$	0.1

3.2.

In vivo Mouse Experiments

In order to further evaluate the performance of the proposed algorithm in practical applications, in vivo mouse experiments were conducted under the protocol approved by the Institutional Animal Care and Use Committee of Tsinghua University. Experiments were performed with the dual-modality FMT/CT system shown in Fig. 1, where the CT images were obtained to examine the reconstruction accuracy. A nude mouse was anesthetized by an abdominal cavity injection of 10% chloral hydrate. In order to mimic the fluorescent target, a transparent glass tube (0.4 cm outer diameter) filled with $1.3 μ M$ ICG was embedded into the abdomen of the mouse. The fusion result of the two dimensional (2-D) fluorescence and white-light images is shown in Fig. 6(a), where the two red dotted lines determine the FMT reconstruction region. The 3-D solid geometry of the mouse is reconstructed from 72 white light images and is shown in Fig. 6(b),¹⁶ where the red line corresponds to the position of the CT transverse reconstruction result shown in Fig. 6(c). In addition, Fig. 6(d) shows the CT-based 3-D rendering of the implanted tube and skeleton inside the mouse surface in order to clarify the morphology of the reconstructed fluorescence tube. Since the tube was embedded in the abdomen of the mouse, the optical properties were assumed to be homogeneous, and the absorption and reduced scattering coefficients were set as $2.38 \times 10^{- 2} {cm}^{- 1}$ and $10.3 {cm}^{- 1}$ , according to Ref. 17.

Fig. 6

In vivo mouse experiments. (a) Fusion of the white-light and fluorescence images. The region between the two red dotted lines is used for FMT reconstruction. (b) 3-D mouse solid geometry reconstructed from 72 white light images, where the red line indicates the position of the reconstructed CT slice (c). (d) Reconstructed CT-based 3-D rendering of the implanted tube and skeleton inside the mouse surface.

In order to compare the quality of the results reconstructed by different L1-norm regularizations with a limited number of projections, the excitation-light and fluorescence-light measurements collected at four projection angles were used in the FMT reconstruction. The dimensions of the corresponding sensitive matrix $W$ were $563 \times 5553$ , and the discretized mesh resolution was set as $0.1 \times 0.1 \times 0.1 {cm}^{3}$ . Figure 7 shows the reconstructed 3-D images and cross sections corresponding to the red line in Fig. 6(b). The first column of Fig. 7 denotes the results reconstructed by the DCT-based reweighted L1-norm regularization and the reconstruction time is 222 seconds. The second column is obtained by the DCT-based L1-norm regularization without the reweighted strategy and the reconstruction time is 274 s. The third and fourth columns correspond to regular L1-norm and classic reweighted L1-norm regularizations, and the corresponding reconstruction times are 122 and 149 s, respectively. Compared with the regular L1-norm and classic L1-norm regularizations, the 3-D DCT transform in the DCT-based regularizations consumes more reconstruction time. The first row of Fig. 7 denotes the fusion FMT/CT slices for different L1-norm regularizations.

Fig. 7

Reconstruction results of the in vivo mouse experiments with four projections. (a1) and (a2) Slice and 3-D results of DCT-based reweighted L1-norm regularization. (b1) and (b2) Slice and 3-D results of DCT-based L1-norm regularization without reweighted strategy. (c1) and (c2) Slice and 3-D results of regular L1-norm regularization. (d1) and (d2) Slice and 3-D results of the classic L1-norm regularization. Each slice corresponds to the position denoted by the red line in Fig. 6(b) and is composed of the reconstructed FMT result fused with the CT slice in Fig. 6(e). Each FMT image is normalized by its maximum.

In order to further demonstrate the localization accuracy of different L1-norm regularizations, the profiles along the white dotted lines in Fig. 7 are shown in Fig. 8. It can be seen that most of the reconstruction-target profiles match the profiles of the true target, though there are some deviations. For this four-projection in vivo FMT reconstruction, the average center-localization error is about 0.09 cm. Deviations of the center of the fluorescent targets reconstructed by different L1-norm regularization strategies may be caused by the inevitable forward-model error and geometrical approximation, such as the rough point-to-point mapping between the captured two dimensional (2-D) CCD images and the light flux on the 3-D mouse surface.

Fig. 8

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

Performances of different L1-norm regularization strategies can be evaluated by the 3-D reconstruction results in Fig. 7. Compared with the result reconstructed by the DCT-based reweighted L1-norm regularization [Fig. 7(a2)], the fluorescent target reconstructed by the DCT-based L1-norm regularization without the reweighted strategy [Fig. 7(b2)] is surrounded by a slight noise on the 3-D surface. This comparison clarifies the fact that the good noise-suppressing ability of the proposed DCT-re-L1 algorithm mainly relies on the proposed DCT-based reweighted strategy rather than the simple DCT-based regularization. Compared with the fluorescent target reconstructed by the DCT-based L1-norm regularization [Figs. 7(a2) and 7(b2)], the target reconstructed by the general L1-norm regularization [Figs. 7(c2) and 7(d2)] has some distortion, which demonstrates the good performance of the DCT-based L1-norm regularization in terms of reconstruction accuracy with regard to the morphology. This is because 3-D DCT can effectively utilize the high degree of space correlation to retain the main morphology of the fluorescent targets. In addition, because DCT can effectively distinguish the noise information, the proposed DCT-based reweighted strategy has a better performance than the classic reweighted strategy in terms of noise suppressing, which can be demonstrated by the comparison between Figs. 7(a2) and 7(d2).

For the in vivo mouse experiments with four projections, the SNR values of the results reconstructed by different L1-norm regularization strategies are listed in Table 3 to quantitatively evaluate the performance of the proposed algorithm. The target and background location was approximated based on the 3-D CT fluorophore distribution, as shown in Fig. 6(d). As shown in Table 3, the result reconstructed by DCT-based reweighted L1-norm regularization has the highest SNR among different L1-norm regularizations, which demonstrates the good performance of the proposed algorithm in terms of noise suppression.

Table 3

SNR of the results reconstructed by different L1-norm regularizations with four projections for the in vivo mouse experiment.

Different L1-norm regularizations	SNR
DCT-based reweighted L1-norm regularization	7.4
DCT-based L1-norm regularization without reweighted strategy	3.9
Regular L1-norm regularization	5.5
Classic reweighted L1-norm regularization	3.6

4. Discussion

In this paper, the FMT inverse problem is formulated as a constrained nonlinear least-square problem with a DCT-based reweighted L1-norm regularization and is solved by the proposed DCT-re-L1 algorithm. In each iteration, the regularization-parameter matrix is adaptively updated based on the values of the DCT coefficients. In order to suppress the noise reconstructed in each iteration, a large regularization parameter (i.e., heavy penalty) is assigned to the small DCT coefficient which usually corresponds to the high-frequency noise. Based on the non-negative region of the previous iterative results, the permission region is adaptively updated, which can expedite the convergence and reduce the memory consumption. In order to clarify the advantages of the proposed DCT-based reweighted L1-norm regularization, different L1-norm regularization methods are employed, including the regular L1-norm regularization, the DCT-based L1-norm regularization without the reweighted strategy and the classic reweighted L1-norm regularization. In order to validate the feasibility of the proposed algorithm for FMT reconstruction, physical phantom and in vivo mouse experiments with a limited number of projections have been conducted. The results demonstrate the good performance of the proposed algorithm in terms of noise suppressing and reconstruction accuracy with regard to the morphology, because 3-D DCT can effectively distinguish the noise information and take advantage of the 3-D space correlation to retain the main morphology of the fluorophores. In addition, in vivo mouse experiments show that the four-projection FMT reconstruction based on the proposed DCT-based reweighted L1-norm regularization can obtain satisfactory reconstruction results with reduced data acquisition time, which contributes to resolving dynamic biological processes in vivo.

As shown in the phantom and in vivo mouse experiments, the 3-D DCT transform causes much time consumption for the FMT reconstruction using the proposed DCT-re-L1 algorithm. In order to further increase the computational speed, fast DCT algorithms (e.g., through fast Fourier transform) can be employed.¹⁸ The in vivo mouse experiments in this paper show that there are some localization errors in the FMT reconstruction. In order to further improve the reconstruction accuracy with regard to localization and morphology for nonconvex fluorescent targets, the forward model with heterogeneous optical properties based on the radiative transfer equation, simplified spherical harmonics approximation or Monte-Carlo is necessary.¹¹^,¹⁹ Additionally, a more accurate mapping model between the 2-D CCD measurements and the 3-D surface photon flux is needed.²⁰ The proposed DCT-based reweighted L1-norm regularization can potentially be utilized in FMT reconstruction with these improved models.

In order to validate the feasibility of the proposed algorithm in practical applications, future work will be focused on in vivo small animal experiments with probe-targeted tumor models. The limited-projection FMT will be applied to the in vivo experiments in order to resolve the tumor metabolism process. In addition, we will incorporate other forms of sparsity-promoting regularization, such as TV, in order to achieve high-quality reconstruction for large (or nonsparse) objects with limited projections.

The focus of this work is the improved SNR using the proposed DCT-re-L1 algorithm when the number of projections is reduced. From the reconstruction results of the phantom experiments, it can be seen that the spatial resolution of the FMT reconstruction declines when the number of projections is reduced. In addition to the number of projections, the reconstruction-image resolution of the FMT system also depends on other experimental parameters such as the arrangement of the excitation-light source, the SNR of the CCD captured images, the thickness of the tissues, the wavelength of the emitted fluorescence light, and the reconstruction algorithms. As pointed out in Ref. 21, compared with L2 regularization, L1 regularization can obtain reconstruction results with higher spatial resolution. The performance of the proposed DCT-re-L1 method in terms of spatial resolution will be investigated 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 Nos. 81227901, 81271617, 61322101, 61361160418 and 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

1.

V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods, 7 603 –614 (2010). http://dx.doi.org/10.1038/nmeth.1483 1548-7091 Google Scholar

2.

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 R1 –R64 (2014). http://dx.doi.org/10.1088/0031-9155/59/1/R1 PHMBA7 0031-9155 Google Scholar

3.

P. Mohajerani et al., “Fluorescence-aided tomographic imaging of synovitis in the human finger,” Radiology, 272 865 –874 (2014). http://dx.doi.org/10.1148/radiol.14132128 RADLAX 0033-8419 Google Scholar

4.

J. Haller et al., “Visualization of pulmonary inflammation using noninvasive fluorescence molecular imaging,” J. Appl. Physiol., 104 795 –802 (2008). http://dx.doi.org/10.1152/japplphysiol.00959.2007 JAPYAA 0021-8987 Google Scholar

5.

J. Shi et al., “Efficient L1 regularization-based reconstruction for fluorescent molecular tomography using restarted nonlinear conjugate gradient,” Opt. Lett., 38 3696 –3699 (2013). http://dx.doi.org/10.1364/OL.38.003696 OPLEDP 0146-9592 Google Scholar

6.

D. Zhu and C. Li, “Nonconvex regularizations in fluorescence molecular tomography for sparsity enhancement,” Phys. Med. Biol., 59 2901 –2912 (2014). http://dx.doi.org/10.1088/0031-9155/59/12/2901 PHMBA7 0031-9155 Google Scholar

7.

A. Behrooz et al., “Total variation regularization for 3D reconstruction in fluorescence tomography: experimental phantom studies,” Appl. Opt., 51 8216 –8227 (2012). http://dx.doi.org/10.1364/AO.51.008216 APOPAI 0003-6935 Google Scholar

8.

D. Wang et al., “High-performance fluorescence molecular tomography through shape-based reconstruction using spherical harmonics parameterization,” PloS One, 9 e94317 (2014). http://dx.doi.org/10.1371/journal.pone.0094317 1932-6203 Google Scholar

9.

W. Xie et al., “Reweighted L1 regularization for restraining artifacts in FMT reconstruction images with limited measurements,” Opt. Lett., 39 4148 –4151 (2014). http://dx.doi.org/10.1364/OL.39.004148 OPLEDP 0146-9592 Google Scholar

10.

N. Ahmed, T. Natarajan and K. R. Rao, “Discrete cosine transform,” IEEE Trans. Comput., C-23 90 –93 (1974). http://dx.doi.org/10.1109/T-C.1974.223784 ITCOB4 0018-9340 Google Scholar

11.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 41 –93 (1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

12.

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 1377 –1386 (2005). http://dx.doi.org/10.1109/TMI.2005.857213 ITMID4 0278-0062 Google Scholar

13.

J. Ripoll et al., “Fast analytical approximation for arbitrary geometries in diffuse optical tomography,” Opt. Lett., 27 527 –529 (2002). http://dx.doi.org/10.1364/OL.27.000527 OPLEDP 0146-9592 Google Scholar

14.

E. J. Candes, M. B. Wakin and S. P. Boyd, “Enhancing sparsity by reweighted L1 minimization,” J. Fourier Anal. Appl., 14 877 –905 (2008). http://dx.doi.org/10.1007/s00041-008-9045-x 1069-5869 Google Scholar

15.

X. Guo et al., “A combined fluorescence and microcomputed tomography system for small animal imaging,” IEEE Trans. Biomed. Eng., 57 2876 –2883 (2010). http://dx.doi.org/10.1109/TBME.2010.2073468 IEBEAX 0018-9294 Google Scholar

16.

D. Wang et al., “In-vivo fluorescence molecular tomography based on optimal small animal surface reconstruction,” Chin. Opt. Lett., 8 82 –85 (2010). http://dx.doi.org/10.3788/COL COLHBT 1671-7694 Google Scholar

17.

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 4225 –4241 (2005). http://dx.doi.org/10.1088/0031-9155/50/17/021 PHMBA7 0031-9155 Google Scholar

18.

E. Feig and S. Winograd, “Fast algorithms for the discrete cosine transform,” IEEE Trans. Signal Process., 40 2174 –21931992 (1992). http://dx.doi.org/10.1109/78.157218 ITPRED 1053-587X Google Scholar

19.

Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express, 17 20178 –20190 (2009). http://dx.doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 Google Scholar

20.

X. Chen et al., “3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images,” Opt. Express, 18 (19), 19876 –19893 (2010). http://dx.doi.org/10.1364/OE.18.019876 OPEXFF 1094-4087 Google Scholar

21.

J. Shi et al., “Enhanced spatial resolution in fluorescence molecular tomography using restarted L1-regularized nonlinear conjugate gradient algorithm,” J. Biomed. Opt., 19 (4), 046018 (2014). http://dx.doi.org/10.1117/1.JBO.19.4.046018 JBOPFO 1083-3668 Google Scholar

Biography

Junwei Shi received his MS degree from Xi’an Jiaotong University, Xi’an, Shaanxi, China, in 2011. Since 2011, he has been working toward his PhD in the Department of Biomedical Engineering, Tsinghua University, Beijing, China. His research interest is fluorescence molecular tomography reconstruction algorithms.

Fei Liu received a bachelor’s degree in biomedical engineering from Zhejiang University, Zhejiang, China, in 2008 and received his PhD in the Department of Biomedical Engineering, Tsinghua University, Beijing, China, in 2013. She is now a postdoctoral research scientist in the School of Medicine, Tsinghua University, Beijing, China. Her research interest is fluorescence molecular tomography for small animal imaging.

Jiulou Zhang is currently a PhD candidate in the Department of Biomedical Engineering at Tsinghua University, Beijing, China. He received a bachelor’s degree in biomedical engineering from Henan University of Science and Technology in 2009 and obtained his master’s degree in biomedical engineering from Southern Medical University in 2012. His research interest is fluorescence molecular tomography for small animal imaging.

Jianwen Luo received his BS and PhD degrees from Tsinghua University in 2000 and 2005, respectively. He was a postdoc and then associate research scientist at Columbia University from 2005 to 2011. He became a professor at Tsinghua University in 2011. He was enrolled in Thousand Young Talents Program of China in 2012 and received Excellent Young Scientists Fund from NSFC in 2013. His research interest includes ultrasound imaging and fluorescence molecular tomography.

Jing Bai received MS and PhD degrees from Drexel University in 1983 and 1985, respectively. From 1985 to 1987, she was a research associate and assistant professor in the Biomedical Engineering and Science Institute, Drexel University. In 1991, she became a professor in the Department of Biomedical Engineering at Tsinghua University. Her current research interests include medical ultrasound and molecular imaging. She has authored or coauthored ten books and more than 300 journal papers.

Citation Download Citation

Junwei Shi, Fei Liu, Jiulou Zhang, Jianwen Luo, and Jing Bai "Fluorescence molecular tomography reconstruction via discrete cosine transform-based regularization," Journal of Biomedical Optics 20(5), 055004 (13 May 2015). https://doi.org/10.1117/1.JBO.20.5.055004

Published: 13 May 2015

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 25 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Reconstruction algorithms

Luminescence

Signal to noise ratio

In vivo imaging

Tomography

3D acquisition

Fluorescence tomography

1.

Introduction

2.

Methods

2.1.

Forward Model for FMT

Eq. (1)

Eq. (2)

2.2.

DCT-Regularized Reconstruction Algorithm

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Algorithm 1

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

3.

Experiments and Results

3.1.

Physical Phantom Experiments

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Table 1

Eq. (15)

Table 2

3.2.

In vivo Mouse Experiments

Fig. 6

Fig. 7

Fig. 8

Table 3

4.

Discussion

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years