Abstract

Through the reconstruction of the fluorescent probe distributions, fluorescence molecular tomography (FMT) can three-dimensionally resolve the molecular processes in small animals in vivo. In this paper, we propose an FMT reconstruction algorithm based on the iterated shrinkage method. By incorporating a surrogate function, the original optimization problem can be decoupled, which enables us to use the general sparsity regularization. Due to the sparsity characteristic of the fluorescent sources, the performance of this method can be greatly enhanced, which leads to a fast reconstruction algorithm. Numerical simulations and physical experiments were conducted. Compared to Newton method with Tikhonov regularization, the iterated shrinkage based algorithm can obtain more accurate results, even with very limited measurement data.

© 2010 OSA

1. Introduction

In vivo small animal optical molecular imaging has become an important and rapidly developing method for biomedical research, and has been widely utilized for cancer detection, drug discovery and gene expression visualization, etc [14]. When small animals are labeled with optical molecular probes, the corresponding biological information can be indirectly visualized by measuring the emitted light photons of the probes over the surface, thus realizing non-invasive investigation of the small animals [5]. Among optical molecular imaging modalities, fluorescence molecular imaging has become a popular and promising technique. This is partly due to the wealth of the fluorescent probes [6]. However, because the photons in the NIR or visible spectrum are highly scattered by tissue, the traditional planar fluorescence imaging cannot reflect the probe distributions accurately. Fluorescence molecular tomography (FMT), on the contrary, can three-dimensionally resolve the molecular processes by measuring the photons over the animal surface and reconstructing the distribution of the fluorescent probes, which makes it a hot spot in recent years.

FMT is often an ill-posed problem since only the photon distribution over the surface is measurable [7]. This can be alleviated by increasing the measurement data sets. However, even if sufficient measurements can be obtained, the problem may still be ill-conditioned, which means that it is unstable and is sensitive to noises caused by CCD measurement errors and data interpolation errors. To compute a meaningful approximate solution, more a priori information should be incorporated to regularize the FMT problem [8]. Among different regularization methods, the Tikhonov regularization is a popular method that has been widely adopted in FMT problems [10,11]. Tikhonov method assumes that the “size” of the solution should not be very large and adds L2-norm constraint of the solution to the original problem, thus making the problem less sensitive to perturbations. The advantage of L2-norm regularization is that the problem is simple and can be solved efficiently by standard minimization algorithms, such as Newton method and conjugate gradient method. However, the solution is often over-smoothed with reduced intensities, and the localized features are lost during the reconstruction process [12]. In recent years, some researchers began to use sparsity regularization for the optical tomography problems [8,9,12]. For FMT this is based on the fact that, the domains of the fluorescent sources are usually very small and sparse compared with the entire reconstruction domain [9]. This can be considered as valuable a priori information for FMT. A straightforward way to integrate the sparsity constraint is to replace the Tikhonov method with L0-norm regularization. However, the optimization problem becomes NP-hard in this case, and cannot be solved efficiently [13]. A tradeoff is to select a proper Lp-norm with p[1,2). When p is within this range, large values in the solution are penalized less severely compared with Tikhonov regularization. Therefore, the Lp regularization (p[1,2)) has a higher tendency to promote the sparsity of the reconstructed result. This effect has been demonstrated in many published articles [8,13]. Another advantage of the sparsity regularization is that, it can still perform well when the measurement data is very limited. This has been well studied in the area of compressed sensing [14]. Besides, the Lp-norm with p[1,2) is convex, so the FMT problem remains convex which is important for successful FMT reconstructions. However, this requires the design of new algorithms that can solve the optimization problems with general Lp-norm regularization. Another problem lies in the performance of the reconstruction algorithm. A long reconstruction time may become an obstacle for FMT to be transferred into practical use. Therefore, we believe that designing a fast reconstruction algorithm is always a hot spot for FMT.

In this paper, a fast reconstruction algorithm based on the iterated shrinkage method is proposed [15]. This algorithm decouples the original high dimensional optimization problem and converts it into a set of one dimensional optimization problems, which enables us to handle each one separately. More importantly, general Lp-norm regularization can be incorporated straightforwardly without increasing the complexity of the problem. Next, we provide two reconstruction strategies for different situations together with their complexity analysis. Last, we explain that duo to the sparsity of the light sources, the computational burden can be greatly reduced, which leads to a fast reconstruction algorithm.

This paper is organized as follows. The iterated shrinkage based reconstruction algorithm is presented in section 2. In section 3, numerical simulations and physical experiments are conducted to evaluate the performance of the proposed method. Finally, we discuss the results and conclude this paper.

2. Method

2.1 Photon propagation model

For photon propagation in biological tissue within the near infrared spectral window, scattering is the dominant phenomenon over absorption. Therefore, the diffuse approximation to the radiative transfer equation has been extensively applied to model the photon transport [1620]. For steady-state FMT with point excitation sources, the following coupled diffuse equations can be utilized to depict the forward problem:

{·(Dx(r)Φx(r))μax(r)Φx(r)=Θδ(rrl)·(Dm(r)Φm(r))μam(r)Φm(r)=Φx(r)ημaf(r)(rΩ)
where subscripts x and m denote the excitation and emission wavelengths, respectively; Ω denotes the domain of the problem; Φx,m is the photon flux density; μax,am is the absorption coefficient and Dx,m=1/3(μax,am+(1g)μsx,sm) is the diffusion coefficient, μsx,sm is the scattering coefficient, and g is the anisotropy parameter; ημaf denotes the fluorescent yield which is to be reconstructed. Here, we assume that the absorption and scattering of the excitation light caused by the fluorescent probes will have little effect on the distribution of Φx, as the fluorescent probes often occupy a very small volume compared with the imaging domain Ω. In this forward model, the excitation light is implemented as isotropic point sources located one mean free path of photon transport beneath the surface at different locationsrl(l=1,2,,L). Θ denotes the amplitude of the point sources. The coupled equations are complemented by Robin-type boundary conditions which depict the refractive index mismatch between the external medium and Ω:
Φx,m(r)+2ADx,m(r)(v(r)·Φx,m(r))(rΩ)
where v denotes the outward normal to the surface. A is a constant depending on the optical reflective index mismatch at the boundary [21].

2.2 Linear relationship establishment

Instead of solving Eqs. (1) and (2) directly, they are posed in the weak solution forms. Discretizing the domain with small elements and employing the base functions as the test functions, the FMT problem can be linearized and the following matrix-form equations can be obtained.

[Kx]{Φx}={Sx}
[Km]{Φm}=[F]{X}
where Kx,m is the system matrix. Matrix F is obtained by discretizing the unknown fluorescent yield distribution. Vector X denotes the fluorescent yield to be reconstructed. For each excitation point source at rl(l=1,2,,L), Φx can be directly obtained by solving Eq. (3). Considering the inverse crime problem, Φx is calculated on a fine mesh using 2nd order Lagrange elements. Then it is projected onto a coarse mesh which will be used for the reconstruction of X with linear elements. As matrix Km is symmetrical positive definite, Eq. (4) can be transformed into:
{Φm,l}=[Km,l1][F]{X}=[Bl]{X}
Removing the unmeasurable entries in Φm,l and corresponding rows in Bl, we have:
{Φm,lmeas}=[Al]{X}
Next, we assemble Eq. (6) for different excitation locations and obtain the following matrix-form equation:
{Φ}=[A]{X}
The detailed descriptions can be found in [7].

2.3 Reconstruction based on iterated shrinkage method

As is mentioned, FMT is usually an ill-posed problem, which means that the dimension of the null space of A is not zero. Therefore, the solution is not unique in this case. Even if the FMT problem becomes less ill-posed when more fluorescence measurement data can be captured by the CCD camera, it can still remain ill-conditioned (with a large condition number). Therefore, errors in the FMT problem can be magnified, which will affect the reconstruction results. Errors are inevitable and can be introduced in several ways, e.g. the fluorescence measurement errors caused by CCD camera noise and the approximation errors caused by data interpolation. A standard way to regularize the problem is to incorporate additional constraints on the solution, which can be considered as a kind of a priori information:

minX0E(X)=12AXΦ22+λR(X)
where R(X) is the penalty function, and λ is a positive real number called the regularization parameter which is used to balance the two terms. In this paper, X is considered to be non-negative. When L2-norm penalty function is used, this becomes the popular Tikhonov regularization method. Here, we only consider the case when R(X) is ||X||pp with p[1,2). In this case, the original energy function for the FMT problem with general Lp-norm regularization can be represented as follows:
E(X)=12AXΦ22+λXpp=12AXΦ22+λi=1Nxip
where N is the dimension of X. Here, we define Xp1=[x1p1,x2p1,,xNp1]T. To minimize the energy function, we calculate the partial derivative of E(X), and set the result to zero:
ATAXATΦ+λpXp1=0
When p=1 and xi=0, Eq. (10) is not correct. This special case will be considered later. Equation (10) is non-linear except when p=1. If ATA is a diagonal matrix, then solving X is equivalent to solving each xi individually, which is quite simple. However, this is generally not the case for FMT, and all xi are closely coupled. Therefore, solving Eq. (10) is not trivial.

Next, we introduce a surrogate function with the following form [15]:

S(X;X0)=c2XX02212AXAX022
where c is a constant scalar and X0 is a constant vector. The parameter c should be chosen so that S(X;X0) is strictly convex. This implies that the Hessian matrix of S(X;X0) should be positive definite, cIATA0. This can be satisfied by choosing parameter c>||ATA||2=ρ(ATA). Adding E(X) and S(X;X0) together, we get the new energy function as follows:
J(X;X0)=12AXΦ22+λXpp+c2XX02212AXAX022
Calculating the partial derivative of Eq. (12) with respect to X and setting the result to zero, we have:
J(X;X0)=[AT(ΦAX0)+cX0]+λpXp1+cX=0
Let D=c1AT(ΦAX0)+X0 represent all the constant items, Eq. (13) can be further simplified as follows:
XD+λpcXp1=0
In Eq. (14), all xi are decoupled now, and solving Eq. (14) is equivalent to solving each xi individually, which makes the use of general Lp-norm regularization possible. This is a distinctive advantage of the iterated shrinkage method compared with other methods, e.g. the Lasso estimator that is based on L1 regularization [22]. Here, we take p=1 for example. To solve a certain xi, we firstly assumes that xi>0. Then, if di>λ/c, we have xi=diλ/c. If di does not satisfy the condition, we don’t calculate xi using Eq. (14) but simply set xi=0. Therefore, we can avoid the aforementioned special case. Now, we define a function named Shrink:
Shrinkλ,c(z;p=1)={zλcz>λc0zλc
This function maps the values smaller than the threshold to zero; values larger than the threshold are “shrinked”, thus the name of the function. An important feature of the Shrink function (not limited to p=1) is that it can be determined once and off-line. Using this Shrink function, the optimal solution to the optimization problem can be represented as:

Xopt=argminXJ(X;X0)=[Shrinkλ,c(d1;p=1)Shrinkλ,c(d2;p=1)Shrinkλ,c(dN;p=1)]

Next, we analyze the influence of the surrogate function on the optimization problem. Since S(X;X0) is strictly convex, it has a unique global minimum. Zeroing the derivative of S(X;X0) with respect to X leads to the equation:

ATA(XX0)=c(XX0)
Since X=X0 is obviously a solution to Eq. (17), it is the unique solution. Therefore, the global minimum of S(X;X0) is zero when X=X0. From the above analysis, we know that the optimal solution to this optimization problem will have a bias towards X0. X0 can be regarded as the a priori knowledge of X. However, this a priori information may not be obtainable in most cases. To resolve this problem, an iterative reconstruction scheme is adopted instead. For a new iteration k+1, Xk+1 is calculated by replacing X0 with Xk:
xk+1,i=Shrinkλ,c(dk,i(Xk);p=1)(i=1,2,,N)
Now, D becomes a function of Xk and can be updated iteratively. Here, we should point out that the convergence analysis in [15] cannot be directly used to analyze the proposed method. A major reason is that, the base functions in the finite element framework are not orthogonal, though they are linearly independent. However, it is shown in reference [23] that even with unorthogonal bases, the iterated shrinkage method is still practiced successfully. We share similar experience in the FMT reconstructions. This can be seen from the reconstruction results in section 3. Here, we provide a brief explanation to show that the proposed method still makes sense in this case. In fact, J(X;Xk) is a upper bound of E(X), which means J(X;Xk)E(X). Equality can be satisfied if and only if X=Xk. Then, for a certain iteration k, the following inequality exists:
E(Xk+1)J(Xk+1;Xk)J(Xk;Xk)=E(Xk)
Therefore, E(Xk) is non-increasing as k grows.

Now, we consider the calculation of c satisfying c>ρ(ATA). A straightforward way is to calculate all the Eigen values of ATA, and c can be set to the largest Eigen value plus ε, where ε is a small positive value. However, this could be rather time-consuming. Therefore, we provide two alternative ways to determine c. For the first one, we estimate the upper bound of the maximum Eigen value instead of actually calculating it, and c can be set to the upper bound plus ε. Several estimation algorithms can be adopted, e.g. the popular Minc method [24]. The second way is to calculate the maximum Eigen value directly using e.g. the Power method [25], and c can be set similarly. Sometimes, the estimated upper bound of maximum Eigen value may be too large compared with the true value. This will increase the influence of the surrogate function on the optimization problem and will slow down the convergence rate. Therefore, for the above two ways, we prefer the latter one, because it is more accurate.

Next, we provide two reconstruction strategies for different situations. Suppose we have obtained A and Φ. For a certain iteration k, Dk can be represented in the following two ways:

Dk=1cAT(ΦAXk)+Xk
Dk=1cATΦ1cATAXk+Xk
Equations (20) and (21) lead to the following two reconstruction strategies:

Tables Icon

Reconstruction strategy 1

For reconstruction strategy 1, we don’t compute ATA since this is time-consuming. As we know, ρ(ATA) is equal to the square of the maximum singular value smax of A, and smax is equal to the maximum Eigen value eigmax of the assembled matrix [0,A;AT,0]. So we can calculate eigmax instead. For each iteration, two matrix-vector multiplications are needed. If we ignore the computations for the vector-vector operations (including the Shrink operations) and the calculation of eigmax, which are relatively cheap in the total reconstruction process, the computational complexity of strategy 1 can be represented as O(kmax*2mn).

For reconstruction strategy 2, we pre-calculate both ATA and ATΦ. The computational complexity of calculating ATA is O(m*n2). For each iteration, only one matrix-vector multiplication is needed. Then, the computational complexity of strategy 2 can be represented as O((kmax+m)*n2).

From the above analysis, it is clear that when kmax is not very large, strategy 1 is more efficient. For a large kmax, especially when kmax is much larger than the dimension of the problem, strategy 2 becomes a better choice.

When kmax is very large, even strategy 2 will be quite time-consuming. Fortunately, if the sparsity characteristic of X is considered, the performance of both strategies can be greatly improved. Suppose Xk is sparse with nk non-zero elements, and nk is usually far less than n. When multiplying Xk by a matrix, the multiplications involving the zero elements in Xk can be simply ignored. Taking strategy 2 for example, the computational complexity of calculating MXk is O(n2). If the sparsity of Xk is considered, only n*nk multiplications are needed, which is a small fraction of the original. For strategy 1, the computational complexity of calculating Vk can be reduced to O(mnk). However, Vk is not sparse any more. Therefore, computing Dk still needs m*n multiplications.

3. Results

3.1 Simulation verifications

In this subsection, heterogeneous simulation experiments were conducted to demonstrate the performance of the proposed method. Figure 1(a) shows the heterogeneous cylindrical phantom we used, which was of 20mm in diameter and 20mm in height. Figure 1(b) is a slice image of the phantom in z = 0 plane. The phantom consisted of four kinds of materials to represent muscle (M), lung (L), heart (H) and bone (B), respectively. The optical parameters can be found in Table 1 . The black dots in Fig. 1(b) represent the excitation light sources, which were modeled as isotropic point sources located one mean free path of photon transport beneath the surface in z = 0 plane.

 

Fig. 1 Mouse-mimicking heterogeneous phantom with four kinds of materials to represent muscle (M), lung (L), heart (H) and bone (B), respectively. (a) is the 3D view of the phantom. (b) is the slice image of the phantom in z = 0 plane. The black dots in (b) represent the excitation point source locations. For each excitation location, fluorescence is measured from the opposite cylindrical side with 160° field of view.

Download Full Size | PPT Slide | PDF

Tables Icon

Table 1. Optical parameters of the heterogeneous phantom [26]

As is mentioned, for FMT, the domains of the fluorescent sources are usually very small and sparse. Therefore, we used small spheres to represent the fluorescent sources in our experiments. Figure 2 . shows the three different phantom setups for single source (left column), double sources (middle column) and three sources (right column), respectively. The first row of Fig. 2 is the 3D views of the phantom setups and the second row is the slice images in z=0 plane. All the fluorescent sources were of 2mm in diameter and were centered in z=0 plane. The fluorescent yields of these sources were set to be 8. These three phantom setups were discretized with tetrahedrons. For the forward problem, three fine meshes, each with about 23000 degrees of freedom (DOFs), were used for different phantom setups. Then the forward solutions were projected onto a single coarse mesh with 2710 DOFs, which was used for the inverse problem. Fluorescence measurement was implemented in transillumination mode. For each excitation source, measurement was taken from the opposite cylindrical side with 160° field of view (FOV), which is illustrated in Fig. 1(b). This means that all the nodes on the cylindrical surface within this FOV were considered to be measurable. When practical fluorescence measurements are taken using a CCD camera, the shot noise always exists, which obeys the Poisson distribution. If large numbers of photons are collected, the shot noise will approach a Gaussian distribution. Therefore, for FMT reconstructions with simulated data, Gaussian noise is often used to simulate the real case. In our simulation experiments, we added 5% Gaussian noise to the measurement data.

 

Fig. 2 Three different phantom setups for single fluorescent source (left column), double fluorescent sources (middle column) and three fluorescent sources (right column), respectively. The first row is the 3D views of the phantoms, and the second row is the slice images in z = 0 plane. All the fluorescent sources are spherical and are centered in z = 0 plane. The diameters of these spherical fluorescent sources are all set to be 2mm.

Download Full Size | PPT Slide | PDF

In this paper, L1-norm regularization was utilized. The maximum iteration number kmax was set to be 30000, which was sufficiently large for these experiments. Since kmax was larger than the dimension of the inverse problem, reconstruction strategy 2 was adopted. To better illustrate the performance of the proposed method, we compared it to the bound-constraint Newton method that minimized the original energy function E(X) with Tikhonov regularization. Since the Newton method is a second-order optimization method, it generally converges faster than those first-order ones. Here, we set the maximum iteration number for the Newton method to be 300. Regularization parameter generally plays an important role in the FMT reconstructions. In this paper, the regularization parameters for the Newton method and the proposed method were manually optimized. Finding the optimal or near-optimal regularization parameters automatically will be our future work. Both reconstruction algorithms were implemented in C++, and all the reconstructions were carried out on a personal computer with 2.33GHz Intel Core2 duo CPU and 2GB RAM.

In the first set of experiments, fluorescence was excited by point sources from 15 different locations in sequence, which is illustrated in Fig. 1(b). Measurements were taken every 24° and a total of 15 data sets were acquired for the reconstruction of the fluorescent yield. Figures 3 , 4 and 5 show the reconstruction results for three different phantom setups, which are presented in the form of slice images in z=0 plane and iso-surfaces for 30% of the maximum value. The small circles in the slice images denote the real positions of the fluorescent sources. From Figs. 3, 4 and 5 we can clearly see that, when only one fluorescent source exists, the results from both methods are quite satisfactory. However, when multiple sources exist, the over-smooth effect of Tikhonov regularization begins to appear. On the contrary, the proposed method can still preserve the sparsity of the sources very well in this case, and the reconstructed intensities are greater. To analyze these results quantitatively, we define the location error to be ||LrL0||2, where L0 is the real location of the source center and Lr is the location of the node with the maximum reconstructed value for that source. We also define the relative intensity error to be |IrI0|/I0, where I0 is the true intensity of the source and Ir is the maximum reconstructed intensity. The quantitative comparisons between these results are presented in Table 2 . From Table 2 we can see that, both reconstruction methods can obtain satisfactory source localizations. But the relative intensity errors for the proposed method are smaller.

 

Fig. 3 Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 1 spherical fluorescent source and 15 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

 

Fig. 4 Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 2 spherical fluorescent sources and 15 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

 

Fig. 5 Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 3 spherical fluorescent sources and 15 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

Tables Icon

Table 2. Quantitative comparisons between the results from the Newton method with Tikhonov regularization (Newton-L2) and the iterated shrinkage based method with L1 regularization (IS-L1) for 15 measurement sets

Next, we compare the running time of the two methods. Since both methods need the pre-calculation of ATA and ATΦ, we set the starting point at the time when ATA and ATΦ were just obtained. Then, the running time of the Newton method for three different phantom setups was 166.62s, 189.32s and 180.09s, respectively. And the running time of the proposed method was 1.16s, 2.11s and 4.46s, respectively. The proposed method is much more efficient.

In the second set of experiments, we reduced the amount of measurement data to simulate a much worse case. This is possible when long-time measurement is not appropriate or feasible. For instance, when imaging small animals like mice, the artifacts caused by movements must be taken into consideration. Besides, long-time measurement can cause the bleaching effect of the fluorescent probe, which may influence the reconstructed results. One way to resolve these problems is to reduce the number of fluorescence measurements. This requires that we should be able to reconstruct the fluorescent sources from very limited data. It has been shown for bioluminescence tomography that, by using sparsity constraint, satisfactory results can still be achievable even with very limited imaging data [8]. Here, we only retained the measurement data sets generated by excitation point source 1, 6 and 11. Figures 6 , 7 and 8 show the reconstruction results in this case. From these figures we can see that, when the measurement data is very limited and multiple fluorescent sources exist, the proposed method can obtain much better results compared to the Newton method with Tikhonov regularization. This demonstrates the applicability of the proposed method under more ill-posed conditions. Quantitative comparisons are presented in Table 3 . For the Newton method, the reconstruction errors for source S1 and S2 in the third phantom setup are not presented, because the two sources cannot be separated in the result. The running time of the Newton method for three different phantom setups was 183.98s, 173.91s and 168.02s, respectively. And the running time of the proposed method was 1.84s, 2.97s and 2.42s, respectively.

 

Fig. 6 Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 1 spherical fluorescent source and 3 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

 

Fig. 7 Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 2 spherical fluorescent sources and 3 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

 

Fig. 8 Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 3 spherical fluorescent sources and 3 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

Tables Icon

Table 3. Quantitative comparisons between the results from the Newton method with Tikhonov regularization (Newton-L2) and the iterated shrinkage based method with L1 regularization (IS-L1) for 3 measurement sets

3.2 Physical experiments

In this subsection, physical experiments were conducted to further test the proposed algorithm. Figure 9 illustrates the experimental setup. Excitation illumination was provided by a 671nm CW laser and the power was set to 20mW. The spot diameter of the laser beam was approximately 1mm. The fluorescence measurements were implemented in transillumination mode. A 10nm band-pass filter centered at 700nm was used to allow light transmission at the emission wave-length. The optical density of the filter at the excitation wave-length was larger than 5. Fluorescence was captured by a CCD camera (Princeton Instruments VersArray 1300B, Roper Scientific, Trenton, NJ) which was cooled to −110°C. To reconstruct the sources, the pixels of the measured fluorescence image should firstly be converted into the corresponding photon flux densities. Then, the calibrated image needs to be back-projected onto the surface of the phantom. Errors can be introduced during the mapping procedure. To avoid this problem and to test only the performance of the proposed reconstruction method, a cubic phantom was utilized in this experiment, and the back-projection was reduced to a point-to-point mapping. Figure 10(a) shows the phantom we used. The side length is 20mm. The phantom is made from polyoxymethylene, and the optical parameters for both excitation and emission wave-lengths, which are presented in Table 4 , were determined by diffuse optical tomography. The sparsity-promoting characteristic of the proposed method can be better demonstrated when more than one source exists, which can be seen from our simulation results. Therefore, in this experiment, two small holes of 1.25mm in radius were drilled to allow two fluorescent sources. 2000nM Cy5.5 solution was used as the fluorescent source. The height of the two cylindrical sources was 2mm, and the centers were at (−3.75mm, 3.75mm, 1mm) and (3.75mm, −3.75mm, 1mm) respectively, which is illustrated in Fig. 10(b). To simulate a badly ill-posed situation, fluorescence was excited by point sources from only 4 different locations in z=0 plane, which is illustrated in Fig. 10(c). The phantom was placed on a rotational stage, which was controlled by computer, to allow measurements from four sides.

 

Fig. 9 The sketch of the experimental setup.

Download Full Size | PPT Slide | PDF

 

Fig. 10 The homogeneous cubic phantom with 2 cylindrical fluorescent sources. (a) is the photograph of the phantom. (b) is the 3D view of the phantom and the sources. (c) is the slice image of the phantom in z = 0 plane. The black dots in (c) represent the excitation point source locations.

Download Full Size | PPT Slide | PDF

Tables Icon

Table 4. Optical parameters of the cubic phantom

In this experiment, the cubic phantom was discretized with tetrahedrons. For the forward problem, a fine mesh with about 21000 DOFs was adopted. And a coarse mesh with 2705 DOFs was used for the inverse problem. Figure 11 shows the reconstruction results which are presented in the form of slice images in z=1mm plane and iso-surfaces for 30% of the maximum value. From Fig. 11 we can clearly see that, due to the badly ill-posed situation and the over-smooth effect of Tikhonov regularization, the reconstructed fluorescent sources from the Newton method are spread within the entire domain, which makes result totally meaningless. The location errors are not presented for the Newton method, because the source locations cannot be identified from the result. On the contrary, the proposed method with L1 regularization can preserve the sparsity of the fluorescent sources very well, though some artifacts exist near the center. The location errors for the fluorescent sources S1 and S2 were 1.21mm and 0.82mm, respectively. The running time of the Newton method and the proposed method was 135.91 seconds and 2.56 seconds, respectively.

 

Fig. 11 Reconstruction results of the cubic phantom from the Newton method (first row) and the iterated shrinkage based method (second row) using 4 measurement data sets. These results are presented in the form of slice images in z = 1mm plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Download Full Size | PPT Slide | PDF

4. Conclusion

In this paper, an iterated shrinkage based reconstruction algorithm is proposed. By integrating an additional surrogate function, the original high dimensional optimization problem can be decoupled into a set of one dimensional ones that can be solved easily. This also enables us to incorporate sparsity regularization in a graceful way. Two reconstruction strategies are provided for different situations together with their complexity analysis. Besides, we explain that due to the sparsity characteristic of the fluorescent sources, the efficiency of the algorithm can be greatly improved, which leads to a fast reconstruction method. Simulation verifications show that the proposed method outperforms the traditional bound-constrained Newton method with Tikhonov regularization in two ways. First, due to the sparsity constraint of our method, obvious improvements can be seen from these reconstruction results. Second, our method is much faster than the Newton method when the fluorescent sources are sparse. This is important if we want to transfer FMT into practical use. We also test our method under a more ill-posed condition and satisfactory results can still be achievable. Reconstruction of experimental data further demonstrates the performance of the proposed method.

For FMT reconstruction, the choice of the regularization parameter will have a significant impact on the results. A large parameter value can make the reconstructed solution deviate from the real distribution, while a small value will have little contribution to the regularization of the problem. Finding the optimal or near-optimal regularization parameter automatically still remains a challenging task. Generally speaking, two strategies can be used: determine the parameter in advance or update it heuristically. This will be our future research.

Another important issue lies in the accuracy of the photon propagation model itself. The diffuse equation has been extensively utilized to describe light transport in biological tissue, yet it is not applicable in certain regions, such as void or more absorptive regions. Several improved models, e.g. higher order approximation to radiative transfer equation, have been proposed to resolve the problem, though more computations are typically needed and the physical meanings are not such explicit. Since FMT reconstruction is a linear inverse problem in nature, the proposed reconstruction algorithm can potentially be utilized in these improved models.

In conclusion, we have developed a fast reconstruction algorithm with sparsity constraint for FMT. Numerical simulations and physical experiments show the merits of our method compared to the Newton method with Tikhonov regularization. In vivo mouse studies using the proposed method will be reported in the future.

Acknowledgments

This paper is supported by the Project for the National Basic Research Program of China (973) under Grant No.2006CB705700, Changjiang Scholars and Innovative Research Team in University (PCSIRT) under Grant No.IRT0645, CAS Hundred Talents Program, CAS scientific research equipment develop program under Grant No. YZ200766, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant No. KGCX2-YW-129, KSCX2-YW-R-262, the National Natural Science Foundation of China under Grant No. 30672690, 30600151, 60532050, 60621001, 30873462, 60910006, 30970769, 30970771, Beijing Natural Science Fund under Grant No.4071003, Science and Technology Key Project of Beijing Municipal Education Commission under Grant No.KZ200910005005.

References and links

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

2. R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001). [PubMed]  

3. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008). [CrossRef]   [PubMed]  

4. J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008). [CrossRef]   [PubMed]  

5. C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002). [CrossRef]   [PubMed]  

6. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006). [CrossRef]   [PubMed]  

7. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=OE-15-26-18300. [CrossRef]   [PubMed]  

8. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source Reconstruction for Spectrally-resolved Bioluminescence Tomography with Sparse a priori Information,” Opt. Express 17(10), 8062–8080 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-10-8062. [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. D. Wang, X. Song, and J. Bai, “Adaptive-mesh-based algorithm for fluorescence molecular tomography using an analytical solution,” Opt. Express 15(15), 9722–9730 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-15-9722. [CrossRef]   [PubMed]  

11. W. Bangerth and A. Joshi, “Adaptive finite element methods for the solution of inverse problems in optical tomography,” Inverse Probl. 24(3), 034011 (2008). [CrossRef]  

12. N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-21-13695. [CrossRef]   [PubMed]  

13. I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm,” IEEE Trans. Signal Process. 45(3), 600–616 (1997). [CrossRef]  

14. E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006). [CrossRef]  

15. I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004). [CrossRef]  

16. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-16-17-13104. [CrossRef]   [PubMed]  

17. Y. Tan and H. Jiang, “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008). [CrossRef]  

18. A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-22-5402. [CrossRef]   [PubMed]  

19. J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-20-15640. [CrossRef]   [PubMed]  

20. C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16(25), 20317–20333 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-25-20317. [CrossRef]   [PubMed]  

21. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef]   [PubMed]  

22. R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,” J. R. Stat. Soc., B 58, 267–288 (1996).

23. M. Elad, B. Matalon, and M. Zibulevsky, “Image denoising with shrinkage and redundant representations,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, New York2, 1924–1931 (2006).

24. H. Minc, Nonnegative matrices, (Wiley, New York, 1988).

25. M. T. Chu and J. L. Watterson, “On a multivariate eigenvalue problem, Part I: Algebraic theory and a power method,” SIAM J. Sci. Comput. 14(5), 1089–1106 (1993). [CrossRef]  

26. A. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-24-9847. [CrossRef]   [PubMed]  

References

  • View by:
  • |
  • |
  • |

  1. 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]
  2. R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001).
    [PubMed]
  3. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
    [CrossRef] [PubMed]
  4. J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
    [CrossRef] [PubMed]
  5. C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002).
    [CrossRef] [PubMed]
  6. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006).
    [CrossRef] [PubMed]
  7. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=OE-15-26-18300 .
    [CrossRef] [PubMed]
  8. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source Reconstruction for Spectrally-resolved Bioluminescence Tomography with Sparse a priori Information,” Opt. Express 17(10), 8062–8080 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-10-8062 .
    [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. D. Wang, X. Song, and J. Bai, “Adaptive-mesh-based algorithm for fluorescence molecular tomography using an analytical solution,” Opt. Express 15(15), 9722–9730 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-15-9722 .
    [CrossRef] [PubMed]
  11. W. Bangerth and A. Joshi, “Adaptive finite element methods for the solution of inverse problems in optical tomography,” Inverse Probl. 24(3), 034011 (2008).
    [CrossRef]
  12. N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-21-13695 .
    [CrossRef] [PubMed]
  13. I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm,” IEEE Trans. Signal Process. 45(3), 600–616 (1997).
    [CrossRef]
  14. E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006).
    [CrossRef]
  15. I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004).
    [CrossRef]
  16. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-16-17-13104 .
    [CrossRef] [PubMed]
  17. Y. Tan and H. Jiang, “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008).
    [CrossRef]
  18. A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-22-5402 .
    [CrossRef] [PubMed]
  19. J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-20-15640 .
    [CrossRef] [PubMed]
  20. C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16(25), 20317–20333 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-25-20317 .
    [CrossRef] [PubMed]
  21. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
    [CrossRef] [PubMed]
  22. R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,” J. R. Stat. Soc., B 58, 267–288 (1996).
  23. M. Elad, B. Matalon, and M. Zibulevsky, “Image denoising with shrinkage and redundant representations,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, New York2, 1924–1931 (2006).
  24. H. Minc, Nonnegative matrices, (Wiley, New York, 1988).
  25. M. T. Chu and J. L. Watterson, “On a multivariate eigenvalue problem, Part I: Algebraic theory and a power method,” SIAM J. Sci. Comput. 14(5), 1089–1106 (1993).
    [CrossRef]
  26. A. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-24-9847 .
    [CrossRef] [PubMed]

2009 (1)

2008 (7)

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
[CrossRef] [PubMed]

J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
[CrossRef] [PubMed]

W. Bangerth and A. Joshi, “Adaptive finite element methods for the solution of inverse problems in optical tomography,” Inverse Probl. 24(3), 034011 (2008).
[CrossRef]

F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-16-17-13104 .
[CrossRef] [PubMed]

Y. Tan and H. Jiang, “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008).
[CrossRef]

J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-20-15640 .
[CrossRef] [PubMed]

C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16(25), 20317–20333 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-25-20317 .
[CrossRef] [PubMed]

2007 (4)

2006 (2)

V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006).
[CrossRef] [PubMed]

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006).
[CrossRef]

2005 (2)

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]

A. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-24-9847 .
[CrossRef] [PubMed]

2004 (2)

I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004).
[CrossRef]

A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-22-5402 .
[CrossRef] [PubMed]

2002 (1)

C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002).
[CrossRef] [PubMed]

2001 (1)

R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001).
[PubMed]

1997 (1)

I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm,” IEEE Trans. Signal Process. 45(3), 600–616 (1997).
[CrossRef]

1996 (1)

R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,” J. R. Stat. Soc., B 58, 267–288 (1996).

1995 (1)

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
[CrossRef] [PubMed]

1993 (1)

M. T. Chu and J. L. Watterson, “On a multivariate eigenvalue problem, Part I: Algebraic theory and a power method,” SIAM J. Sci. Comput. 14(5), 1089–1106 (1993).
[CrossRef]

Adibi, A.

Arridge, S. R.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
[CrossRef] [PubMed]

Bachmann, M. H.

C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002).
[CrossRef] [PubMed]

Bai, J.

Bangerth, W.

Bao, S.

J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
[CrossRef] [PubMed]

Candès, E. J.

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006).
[CrossRef]

Cao, N.

Chan, T. F.

Chatziioannou, A. F.

Chen, N.

Chu, M. T.

M. T. Chu and J. L. Watterson, “On a multivariate eigenvalue problem, Part I: Algebraic theory and a power method,” SIAM J. Sci. Comput. 14(5), 1089–1106 (1993).
[CrossRef]

Cong, A. X.

Contag, C. H.

C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002).
[CrossRef] [PubMed]

Daubechies, I.

I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004).
[CrossRef]

Defrise, M.

I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004).
[CrossRef]

Delpy, D. T.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
[CrossRef] [PubMed]

DeMol, C.

I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004).
[CrossRef]

Dinkelborg, L. M.

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
[CrossRef] [PubMed]

Douraghy, A.

Eftekhar, A. A.

Feng, J.

Gambhir, S. S.

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
[CrossRef] [PubMed]

Gao, F.

Gorodnitsky, I. F.

I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm,” IEEE Trans. Signal Process. 45(3), 600–616 (1997).
[CrossRef]

Hiraoka, M.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
[CrossRef] [PubMed]

Huang, J.

Jacobs, M.

Jia, K.

Jiang, H.

Y. Tan and H. Jiang, “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008).
[CrossRef]

Joshi, A.

Li, Y.

J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
[CrossRef] [PubMed]

Liang, W.

J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
[CrossRef] [PubMed]

Liu, K.

Lu, Y.

Lv, Y.

Mahmood, U.

R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001).
[PubMed]

Marjono, A.

Mohajerani, P.

Nehorai, A.

Ntziachristos, V.

V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006).
[CrossRef] [PubMed]

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]

Qin, C.

Rao, B. D.

I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm,” IEEE Trans. Signal Process. 45(3), 600–616 (1997).
[CrossRef]

Ripoll, J.

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]

Romberg, J. K.

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006).
[CrossRef]

Schweiger, M.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
[CrossRef] [PubMed]

Sevick-Muraca, E. M.

Song, X.

Stout, D.

Tan, Y.

Y. Tan and H. Jiang, “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008).
[CrossRef]

Tanikawa, Y.

Tao, T.

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006).
[CrossRef]

Tian, J.

Tibshirani, R.

R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,” J. R. Stat. Soc., B 58, 267–288 (1996).

van Bruggen, N.

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
[CrossRef] [PubMed]

Wang, D.

Wang, G.

Wang, H.

Wang, L. V.

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]

Watterson, J. L.

M. T. Chu and J. L. Watterson, “On a multivariate eigenvalue problem, Part I: Algebraic theory and a power method,” SIAM J. Sci. Comput. 14(5), 1089–1106 (1993).
[CrossRef]

Weissleder, R.

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]

R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001).
[PubMed]

Willmann, J. K.

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
[CrossRef] [PubMed]

Xu, M.

Yamada, Y.

Yan, G.

Yan, X. P.

J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
[CrossRef] [PubMed]

Yang, X.

Zhang, L.

Zhang, X.

Zhao, H.

Zhu, S.

Annu. Rev. Biomed. Eng. (2)

C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002).
[CrossRef] [PubMed]

V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006).
[CrossRef] [PubMed]

Appl. Opt. (1)

Commun. Pure Appl. Math. (2)

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006).
[CrossRef]

I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004).
[CrossRef]

IEEE Eng. Med. Biol. Mag. (1)

J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008).
[CrossRef] [PubMed]

IEEE Trans. Signal Process. (1)

I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm,” IEEE Trans. Signal Process. 45(3), 600–616 (1997).
[CrossRef]

Inverse Probl. (1)

W. Bangerth and A. Joshi, “Adaptive finite element methods for the solution of inverse problems in optical tomography,” Inverse Probl. 24(3), 034011 (2008).
[CrossRef]

J. R. Stat. Soc., B (1)

R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,” J. R. Stat. Soc., B 58, 267–288 (1996).

Med. Phys. (2)

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
[CrossRef] [PubMed]

Y. Tan and H. Jiang, “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008).
[CrossRef]

Nat. Biotechnol. (1)

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]

Nat. Rev. Drug Discov. (1)

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008).
[CrossRef] [PubMed]

Opt. Express (9)

X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=OE-15-26-18300 .
[CrossRef] [PubMed]

Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source Reconstruction for Spectrally-resolved Bioluminescence Tomography with Sparse a priori Information,” Opt. Express 17(10), 8062–8080 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-10-8062 .
[CrossRef] [PubMed]

A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-22-5402 .
[CrossRef] [PubMed]

J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-20-15640 .
[CrossRef] [PubMed]

C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16(25), 20317–20333 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-25-20317 .
[CrossRef] [PubMed]

N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-21-13695 .
[CrossRef] [PubMed]

D. Wang, X. Song, and J. Bai, “Adaptive-mesh-based algorithm for fluorescence molecular tomography using an analytical solution,” Opt. Express 15(15), 9722–9730 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-15-9722 .
[CrossRef] [PubMed]

F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-16-17-13104 .
[CrossRef] [PubMed]

A. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-24-9847 .
[CrossRef] [PubMed]

Radiology (1)

R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001).
[PubMed]

SIAM J. Sci. Comput. (1)

M. T. Chu and J. L. Watterson, “On a multivariate eigenvalue problem, Part I: Algebraic theory and a power method,” SIAM J. Sci. Comput. 14(5), 1089–1106 (1993).
[CrossRef]

Other (2)

M. Elad, B. Matalon, and M. Zibulevsky, “Image denoising with shrinkage and redundant representations,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, New York2, 1924–1931 (2006).

H. Minc, Nonnegative matrices, (Wiley, New York, 1988).

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)

Fig. 1
Fig. 1

Mouse-mimicking heterogeneous phantom with four kinds of materials to represent muscle (M), lung (L), heart (H) and bone (B), respectively. (a) is the 3D view of the phantom. (b) is the slice image of the phantom in z = 0 plane. The black dots in (b) represent the excitation point source locations. For each excitation location, fluorescence is measured from the opposite cylindrical side with 160° field of view.

Fig. 2
Fig. 2

Three different phantom setups for single fluorescent source (left column), double fluorescent sources (middle column) and three fluorescent sources (right column), respectively. The first row is the 3D views of the phantoms, and the second row is the slice images in z = 0 plane. All the fluorescent sources are spherical and are centered in z = 0 plane. The diameters of these spherical fluorescent sources are all set to be 2mm.

Fig. 3
Fig. 3

Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 1 spherical fluorescent source and 15 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Fig. 4
Fig. 4

Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 2 spherical fluorescent sources and 15 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Fig. 5
Fig. 5

Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 3 spherical fluorescent sources and 15 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Fig. 6
Fig. 6

Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 1 spherical fluorescent source and 3 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Fig. 7
Fig. 7

Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 2 spherical fluorescent sources and 3 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Fig. 8
Fig. 8

Reconstruction results from the Newton method (first row) and the iterated shrinkage based method (second row) for 3 spherical fluorescent sources and 3 measurement data sets. These results are presented in the form of slice images in z = 0 plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Fig. 9
Fig. 9

The sketch of the experimental setup.

Fig. 10
Fig. 10

The homogeneous cubic phantom with 2 cylindrical fluorescent sources. (a) is the photograph of the phantom. (b) is the 3D view of the phantom and the sources. (c) is the slice image of the phantom in z = 0 plane. The black dots in (c) represent the excitation point source locations.

Fig. 11
Fig. 11

Reconstruction results of the cubic phantom from the Newton method (first row) and the iterated shrinkage based method (second row) using 4 measurement data sets. These results are presented in the form of slice images in z = 1mm plane (left column) and iso-surfaces for 30% of the maximum value (right column). The small circles in the slice images denote the real positions of the fluorescent sources.

Tables (5)

Tables Icon

Table 1 Reconstruction strategy 1

Tables Icon

Table 1 Optical parameters of the heterogeneous phantom [26]

Tables Icon

Table 2 Quantitative comparisons between the results from the Newton method with Tikhonov regularization (Newton-L2) and the iterated shrinkage based method with L1 regularization (IS-L1) for 15 measurement sets

Tables Icon

Table 3 Quantitative comparisons between the results from the Newton method with Tikhonov regularization (Newton-L2) and the iterated shrinkage based method with L1 regularization (IS-L1) for 3 measurement sets

Tables Icon

Table 4 Optical parameters of the cubic phantom

Equations (21)

Equations on this page are rendered with MathJax. Learn more.

{ · ( D x ( r ) Φ x ( r ) ) μ a x ( r ) Φ x ( r ) = Θ δ ( r r l ) · ( D m ( r ) Φ m ( r ) ) μ a m ( r ) Φ m ( r ) = Φ x ( r ) η μ a f ( r ) ( r Ω )
Φ x , m ( r ) + 2 A D x , m ( r ) ( v ( r ) · Φ x , m ( r ) ) ( r Ω )
[ K x ] { Φ x } = { S x }
[ K m ] { Φ m } = [ F ] { X }
{ Φ m , l } = [ K m , l 1 ] [ F ] { X } = [ B l ] { X }
{ Φ m , l m e a s } = [ A l ] { X }
{ Φ } = [ A ] { X }
min X 0 E ( X ) = 1 2 A X Φ 2 2 + λ R ( X )
E ( X ) = 1 2 A X Φ 2 2 + λ X p p = 1 2 A X Φ 2 2 + λ i = 1 N x i p
A T A X A T Φ + λ p X p 1 = 0
S ( X ; X 0 ) = c 2 X X 0 2 2 1 2 A X A X 0 2 2
J ( X ; X 0 ) = 1 2 A X Φ 2 2 + λ X p p + c 2 X X 0 2 2 1 2 A X A X 0 2 2
J ( X ; X 0 ) = [ A T ( Φ A X 0 ) + c X 0 ] + λ p X p 1 + c X = 0
X D + λ p c X p 1 = 0
S h r i n k λ , c ( z ; p = 1 ) = { z λ c z > λ c 0 z λ c
X o p t = arg min X J ( X ; X 0 ) = [ S h r i n k λ , c ( d 1 ; p = 1 ) S h r i n k λ , c ( d 2 ; p = 1 ) S h r i n k λ , c ( d N ; p = 1 ) ]
A T A ( X X 0 ) = c ( X X 0 )
x k + 1 , i = S h r i n k λ , c ( d k , i ( X k ) ; p = 1 ) ( i = 1 , 2 , , N )
E ( X k + 1 ) J ( X k + 1 ; X k ) J ( X k ; X k ) = E ( X k )
D k = 1 c A T ( Φ A X k ) + X k
D k = 1 c A T Φ 1 c A T A X k + X k

Metrics