Open Access
21 April 2017 Laplacian manifold regularization method for fluorescence molecular tomography
Xuelei He, Xiaodong Wang, Huangjian Yi, Yanrong Chen, Xu Zhang, Jingjing Yu, Xiaowei He
Author Affiliations +
Abstract
Sparse regularization methods have been widely used in fluorescence molecular tomography (FMT) for stable three-dimensional reconstruction. Generally, 1-regularization-based methods allow for utilizing the sparsity nature of the target distribution. However, in addition to sparsity, the spatial structure information should be exploited as well. A joint 1 and Laplacian manifold regularization model is proposed to improve the reconstruction performance, and two algorithms (with and without Barzilai–Borwein strategy) are presented to solve the regularization model. Numerical studies and in vivo experiment demonstrate that the proposed Gradient projection-resolved Laplacian manifold regularization method for the joint model performed better than the comparative algorithm for 1 minimization method in both spatial aggregation and location accuracy.

1.

Introduction

Fluorescence molecular tomography (FMT)1,2 is a noninvasive molecular imaging technique. In cooperation with appropriate fluorescent probes, FMT allows for three-dimensional (3-D) imaging of molecular functions and the visualization of biological and physiological processes, which make it powerful in clinical and preclinical research, e.g., cancer and drug development.3,4 To locate and quantify fluorescence probes from the measured light intensities on the surface of biological tissues, both accurate transportation model and sophisticated inverse algorithm are indispensable. However, the inherently ill-posed nature of the inverse problem along with complex scattering and absorbing in the process of near-infrared photons propagation in biological tissues makes FMT challenging.5,6

To generate a meaningful and stable solution, many regularization methods have been introduced to solve the inverse problem. In the early stage of FMT development, 2-norm regularization method, such as Tikhonov regularization,7,8 was widely used to deal with the ill-posedness of FMT. However, 2-norm penalty often leads to blurred reconstructed images due to its over-smoothness. To improve the imaging quality of FMT, various kinds of sparse regularization methods with sparsity-inducing norms have been applied. Theoretically, 0-norm model is expected to produce the sparsest solution, but it is computationally infeasible for practical application. Although greedy algorithms can approximately solve 0-norm model, they are easily trapped into local optimal solution.9,10 As an alternative solution, 1-norm regularization is the most popular due to the existence of plenty of efficient algorithms to solve it.1114 Recently, nonconvex p-norm regularization methods have been investigated in the field of FMT. Compared with 1-norm methods, nonconvex p-norm methods can provide more accurate location and sparser solution.1518

Generally, regularization items are introduced into inverse problem for special goals. For example, considering that the spatial distribution of biological activity studied by FMT tends to occur in localized regions, sparsity-inducing p-norm (p[0,1]) regularizers are utilized to deal with the ill-posedness of inverse problem. Similarly, total variation (TV) penalty is used to promote smoothness of solution while preserving edges of targets.1922 Although TV regularization has been successfully used in image processing, there are some obstacles in applying TV-norm to FMT. As the finite element mesh used for tomography is not as regular as image pixels, it brings heavy computation and great memory cost to optimize the corresponding TV models.20 Except for TV regularization, Laplacian regularization has been used to incorporate structural information of x-ray computed tomography\magnetic resonance to the inverse problem in FMT modalities.2325

Based on the assumption that observed biological mechanisms or activities in FMT are confined to small regions, 1-norm regularization-based methods are widely used in FMT reconstruction. Nevertheless, the 3-D spatial distribution of the fluorescence sources also relies on the tracers or fluorescence probes used in FMT. The source distributions are expected to be localized for highly specific tracers, whereas the spatial structure of nonspecific tracers such as indocyanine green may be more globally distributed.26 For such cases, a penalty term taking spatial structure information into account can be expected to outperform 1-type penalties. Consequently, an inverse model combining both sparsity and spatial structure information is potentially beneficial for recovery of fluorescent targets with comprehensive consideration of the changing situations in FMT applications.

To utilize the spatial structure information, manifold-based learning method is needed. Since Isomap27 and locally linear embedding28 were proposed, manifold learning has become one of the hottest research fields in machine learning. It has been widely applied in semisupervised classification,29 clustering,30 dimensionality reduction,31,32 and so on. Inspired by manifold learning, we introduce Laplacian manifold regularization to obtain spatial aggregation and combine it with sparse regularization to improve FMT reconstruction. Moreover, to solve the sparse reconstruction model with Laplacian manifold regularization, we present two gradient projection algorithms (with and without Barzilai–Borwein strategy), gradient projection-resolved Laplacian manifold (GPRLM)-basic and GPRLM-BB, which can be regarded as the extension of gradient projection for sparse reconstruction (GPSR)33 method for 1-norm regularization model.

The originality of this work lies in the application of Laplacian manifold regularization to FMT and the efficient algorithms to deal with the joint 1-norm and Laplacian manifold regularization model. To assess the performance of the proposed inverse model and reconstruction algorithms, we conducted numerical simulations and in vivo experiment to compare the proposed algorithms with the corresponding version of GPSR.

Rest of this paper is arranged as follows: in Sec. 2, we briefly review the photon propagation model, explain the joint 1-norm and Laplacian manifold regularization model in detail, and present the two algorithms. In Sec. 3, we assess the proposed model and algorithms with simulations and experiments. In Sec. 4, we conclude the paper with a discussion of the results.

2.

Methods

2.1.

Photon Propagation Model

In this paper, we use diffusion equation (DE) for modeling photon propagation in excitation and emission process.34,35 The coupled DEs are given as

Eq. (1)

{·(Dxϕx)μaxϕx=Θsδ(rrsl),·(Dmϕm)μamϕm=ϕxημaf,
where r is the position vector; D,ϕ,μ, and η are the function of r. ϕx(r) and ϕm(r) denote photon intensity (photons/cm2s) at position r, respectively, in which the subscripts x and m represent the excitation and emission processes. Dx,m is the optical diffusion coefficient, which equals to 1/3(μax,m+(1g)μsx,m), in which g represents the anisotropy parameter, μax,m is the absorption coefficient, and μsx,m is the scattering coefficient. η(r) denotes the fluorescent yield efficiency and μaf is the absorption coefficient of fluorophores. rsl is the position of different excitation sources with the amplitude Θs. The boundary conditions accompanying Eq. (1) are modified Robin-type boundary conditions expressed as

Eq. (2)

{n·(Dxϕx)+αxϕx=0,n·(Dmϕm)+αmϕm=0,
where n is the normal vector of the boundary. αx and αm are the Robin boundary coefficients.

Using the finite element method (FEM),36,37 two matrix equations are yielded as

Eq. (3)

{Kxϕx=Lx,Kmϕm=FX,
where Eq. (3) corresponds to Eq. (1), respectively. For a single excitation source with the position rsl and the amplitude Θs, we can calculate the ϕx,sl with Eq. (3) and then ϕm,sl with Eq. (3), namely

Eq. (4)

{ϕx,sl=Kx,sl1Lxϕm,sl=Km,sl1FX=AslX,
where Asl=Km,sl1F.

After removing the measurements that cannot be observed, we composite all systems of linear equations in Eq. (4). The final matrix equation is formed as

Eq. (5)

AX=ϕ,
where X, an N×1 vector, represents distribution of the unknown fluorescent source; ϕ, an M×1 vector, is the boundary observation, and A denotes an M×N weight matrix depending on the geometry and optical parameters.

2.2.

Joint 1 and Laplacian Manifold Regularization Model

The aim of FMT can be boiled down to finding the distribution of the unknown fluorescent target X by solving Eq. (5). However, it is inaccurate to solve the equation directly due to the high ill-posed nature and the existence of noise. To obtain stable and accurate solution of Eq. (5), some form of regularization method is indispensable in the inversion. The widely used 1-norm regularization is to formulate the inverse problem into the following optimization with 1-norm penalty:

Eq. (6)

argminX12ϕAX22+τX1.

Unlike the previous reconstruction method for FMT, we introduce manifold regularization into the inversion process to further improve the quality of reconstruction. Specifically, we convert Eq. (5) into an optimization problem with joint 1 and Laplacian manifold terms

Eq. (7)

argminX12ϕAX22+τX1+λ2XTLX,
where XTLX is the Laplacian manifold regularization term and L is the graph Laplacian matrix, which will be discussed in detail below. Equation (7) is constructed based on a graph model, which is derived by the finite element mesh used for reconstruction. Figure 1 shows the relationship between mesh and graph model, where a 2-D finite element mesh is used for illustrative purposes only. The definition of the graph model is as follows.

Fig. 1

Illustration of graph model based on 2-D finite element mesh, where the nodes and edges of the mesh serve as the nodes and edges of the graph model, respectively.

JBO_22_4_045009_f001.png

The nodes and edges of FEM mesh serve as the vertexes and edges of the graph model, respectively. Let X=(x1,x2,,xN)T be the distribution of fluorescent intensity and pi(i=1,2,,N) be the space coordinate vector of the i’th node corresponding to xi. Then, eij=1 (i,j=1,2,,N) represents that there is an edge between pi and pj (see the red line segment and the red nodes in Fig. 1). Otherwise, if there is no edge between pk and pl (see the blue nodes in Fig. 1), ekl=0 (k,l=1,2,,N).

The Laplacian manifold regularizer is defined as

Eq. (8)

i=1Nj=1Nwij(xixj)2,
where wij is the affine weight of the edge between node pi and node pj. To measure the connection of different points, the weight wij is defined as

Eq. (9)

wij={exp(pipj22σ2)if  eij=1,0if  i=joreij=0.,
where σ>0 is the parameter to adjust the weight matrix. It is apparent that wij[0,1], and the closer pi from pj, the closer wij to 1. By simple mathematical derivation,32 we can obtain a plain form of Eq. (7)

Eq. (10)

i=1Nj=1Nwij(xixj)2=2XT(DW)X,
where W is a symmetrical weight matrix, which equivalents to (wij)N×N. D is a diagonal matrix, which equals to diag(d1,d2,,dN), where di=j=1Nwij=k=1Nwki. Let L=DW, then we obtain the Laplacian manifold regularization item of Eq. (7).

From the definition of Laplacian manifold regularization item of Eq. (7), it is easy to find that (xixj)2 will contribute much to the objective function Eq. (7) if there is a relatively large difference between xi and xj. Therefore, the Laplacian manifold regularizer encourages the adjacent nodes (eij=1) to have similar intensity value and hence it induces the solution to present spatial aggregation. To solve Eq. (7), we extended the GPSR33 method to fit the new joint regularization model.

2.3.

Gradient Projection Method for Joint Regularization

In a similar way as what was done by Figueiredo et al.,33 the variable X in Eq. (7) can be split into its positive and negative parts Zf=(xi)+ and Zl=(xi)+ for all i=1,2,,N, where (·)+ denotes the “positive-part operator.” Then, the Eq. (7) can be reformulated as following bound-constrained quadratic program:

Eq. (11)

argminZf,Zl12ϕA(ZfZl)22+τ1NTZf+τ1NTZl+λ2(ZfZl)TL(ZfZl),s.t.  Zf0,Zl0,
where 1N is length-N vector with 1 as its elements. By substitution, an equivalent quadratic program with non-negative constraints can be obtained

Eq. (12)

argminZ12ZTBZ+CTZ,s.t.  Z0,
where
Z=[ZfZl],C=[τ1NATϕτ1N+ATϕ]
and

Eq. (13)

B=[ATA+λL(ATA+λL)(ATA+λL)ATA+λL].

Following the studies, we present two algorithms to solve the quadratic problem in Eq. (12) based on basic gradient projection method and the Barzilai–Borwein gradient projection method. The two algorithms are denoted as GPRLM-basic and GPRLM-BB, respectively.

For simplicity, let F(Z)=12ZTBZ+CTZ. With gradient projection technique, the solution of Eq. (12) evolves from Z(k) to Z(k+1) in the following way:

Eq. (14)

Zk+1=Z(k)+ξk[W(k)Z(k)],
where W(k)=[Z(k)α(k)F(Z(k))]+.

The difference between GPRLM-basic and GPRLM-BB lies in their searching direction and searching step size. More specifically, their choices of the parameters α(k) and ξ(k). In the basic approach GPRLM-basic, the searching is along the negative gradient F(Z(k)) and confined in the non-negative orthant. Backtracking line search is performed to ensure F(Z) decreases at every iteration. The searching direction for GPRLM-BB algorithm is based on Newton method with approximate Hessian matrixes. It is different from the basic algorithm, which is defined by steepest descent strategy. GPRLM-basic algorithm uses inexact line searching to find the searching step size α(k) and GPRLM-BB algorithm calculates the step size by ξ(k)=mid(0,(δ(k))TF(Z(k))γ(k),1), which can ensure the convergence of GPRLM-BB algorithm.

We define the vector g(k) as

Eq. (15)

gi(k)={[F(Z(k))]i,if  zi(k)>0or[F(Z(k))]i<0,0,otherwise,
where Z(k) is the k’th iteration point obtained by GPRLM-basic, and [F(Z(k))]i is the i’th element of F(Z(k)). Then, F(Z)=BZ+C. The scalar parameter, α, used to yield the step along the gradient descent direction. The initial value of α(k) is chose to be α0=argminαF(Z(k)αg(k)), which can be computed as

Eq. (16)

α0=(g(k))Tg(k)(g(k))TBg(k).

In addition, the operator mid(a,b,c) is defined to return the middle value of {a,b,c}, which is used to confine parameter α(k) and ξ(k) to a given interval.

The detailed flow of algorithms is presented in Fig. 2.

Fig. 2

Flow diagram of GPRLMs. (a) GPRLM-basic and (b) GPRLM-BB.

JBO_22_4_045009_f002.png

3.

Experiments

To assess the joint regularization model and the proposed algorithms, we conducted simulations on a 3-D digital mouse model and in vivo mouse experiment. Since the focus of this section is to validate the effectiveness of Laplacian regularization, the original algorithms GPSR-basic and GPSR-BB proposed by Figueiredo et al.33 instead of more sophisticated algorithms are chosen as competitive algorithms. Consequently, there are four algorithms are investigated in the following simulations and experiments.

For the following simulations and experiments, we used the measurements on all of the boundary mesh nodes for reconstruction. For a fair comparison, we used the same regularization parameter τ for 1-norm penalty term used in both standalone 1 regularization model and the joint regularization model. In addition, the regularization parameter τ was automatically selected by L-curve38 and the parameter λ was selected manually. Specifically, parameter τ ranged from 1e07 to 1e10 and parameter λ for Laplacian regularization term ranged from 1e01 to 1e03. The iteration numbers for all of the testing algorithms were set to 500. We set the 5% of maximum observed fluorescent intensity as the noise standard deviation in simulation.

3.1.

Evaluation Criterion

To quantitatively assess reconstruction, four metrics are used to compare the results, including location error (LE), contrast-to-noise ratio (CNR), and reconstructed time (time). LE is the Euclidean distance between the centers of the reconstructed target and the actual target, which is a measure of location accuracy. CNR39 is defined as the difference between the averaged optical coefficient within the region of interest and the difference within the background region, divided by averaged optical coefficient variation in the background. It is formulated as

Eq. (17)

CNR=μnoiμnob(wnoiσnoi2+wnobσnob2)1/2,
where μnoi is the mean value of the fluorescent yield of the interest nodes and μnob is the mean value of the fluorescent yield of the background nodes. wnoi and wnob are the percentages of the interest nodes and the background nodes in the whole area, respectively. σnoi2 and σnob2 are the standard variance of the interest nodes and the background nodes, respectively.

3.2.

Numerical Simulation on a Three-Dimensional Digital Mouse Model

We conducted several groups of numerical simulations on a 3-D digital mouse model to test the proposed algorithms in terms of the above metrics. The 3-D digital mouse atlas40 is used to offer the structural information and only 340 slices of the mouse atlas are chosen to be investigated, which is composed by six organs including muscle, heart, liver, stomach, kidneys, and lungs. The optical parameters of organs are shown in Table 1.41

Table 1

Optical parameters of the mouse organs.

Organ650 nm670 nmg
μax(mm−1)μsx′(mm−1)μam(mm−1)μsm′(mm−1)
Muscle0.00521.080.00681.030.9
Heart0.00831.010.01040.990.85
Liver0.03290.700.01760.650.9
Stomach0.01141.740.00701.360.92
Kidneys0.00662.250.03802.200.86
Lungs0.01331.970.02031.950.9
Fluorophore0.32921.00.03814.740.9

3.2.1.

Reconstruction of single target

The first group of simulation is to reconstruct a single cylindrical target with 2-mm diameter and 2-mm height located in the liver of the mouse model, with the center at (12,11.5,16) mm. Figure 3(a) shows the mouse model and the target setting. Figure 3(b) shows the distribution of 18 excitation sources located in the plane of Z=16  mm, where the dots denote positions of the excitation sources. The measurement acquisition was performed at a field of view (FOV) of 120 deg opposite the excitation source.

Fig. 3

(a) The mouse model with a cylinder fluorescent target in liver, (b) the positions of 18 excitation sources in the slice of Z=16  mm and the illustrations of a measurable 120 deg FOV on surface opposite excitation sources, and (c) the forward mesh and the simulated distribution of photons density on surface.

JBO_22_4_045009_f003.png

To obtain simulated measurement, forward calculation was carried out on a forward mesh with 31,648 nodes and 177,425 tetrahedrons. Using FEM solving the forward model with known optical absorption and scattering parameters, we can obtain the simulated boundary measurement, as shown in Fig. 3(c). 3-D reconstructions were performed on a tetrahedral mesh discretizing the mouse model with 3020 nodes and 14,903 tetrahedrons.

Figure 4 shows the cross-sectional views of the reconstruction results at Z=16  mm by GPSRs and GPRLMs. Table 2 lists the detailed single-target results by the four comparative algorithms. The reconstruction quality by GPRLMs is better than GPSRs in intensity by visual inspection. The proposed GPRLMs performed slightly better in terms of LE and CNR, as shown in Fig. 4 and Table 2.

Fig. 4

Reconstruction results of single-target experiments. (a), (b), (c), and (d) are reconstructed targets in slice Z=16  mm obtained by GPSR-basic, GPRLM-basic, GPSR-BB, and GPRLM-BB, respectively, (e)–(h) are reconstructed targets in 3-D views corresponding to (a)–(d), respectively.

JBO_22_4_045009_f004.png

Table 2

Comparison results in single-target simulations.

MethodReconstruction center (mm)CNRLE (mm)Time (s)
GPSR-basic(12.25, 11.61, 16.43)17.080.5176.35
GPRLM-basic(12.22, 11.59, 16.44)17.550.5077.02
GPSR-BB(12.14, 11.48, 16.47)16.400.5139.79
GPRLM-BB(12.16, 11.50, 16.47)16.700.5044.13
Note: The best results are in bold.

3.2.2.

Assessment of algorithm stability

Stability analysis was also performed to test the influence of several factors on the proposed reconstruction algorithms with single-target simulation. We considered the position of excitation sources plane, the number of excitation sources, and different levels of noise, respectively. The mouse model and simulation settings were the same as in Sec. 3.2.1. To test the influence of position of excitation sources plane, we moved excitation plane along the Z-axis from the plane Z0 to Z2 (see Fig. 5). In addition, we kept the excitation plane fixed at Z=16  mm and gradually reduced the number of excitation sources or added different levels of noise to observe the change of reconstruction quality.

Fig. 5

Simulation settings: (a) and (b) are XZ and YZ views of different excitation planes, where Z0=16  mm, Z1=12  mm, and Z2=20  mm.

JBO_22_4_045009_f005.png

Figure 6 shows the quantitative results by the comparative algorithms in terms of CNR, LE, and time under different excitation planes. GPRLMs performed better than GPSRs in spatial aggregation in all three test cases from Fig. 6(a). It can be observed that GPRLMs possessed the best location accuracy for all of the tested excitation positions in Fig. 6(b). In Fig. 6(c), it is obvious that BB algorithms (GPSR-BB and GPRLM-BB) run faster than basic algorithms (GPSR-basic and GPRLM-basic).

Fig. 6

Comparison results of different excitation node positions: (a) CNR, (b) LE, and (c) time.

JBO_22_4_045009_f006.png

Figure 7 shows the comparison results with a different number of excitation sources by different reconstruction algorithms. As expected, the location accuracy of the reconstructions and the reconstruction time did witness significant deterioration with the decrease of excitation nodes from Fig. 7(b). Generally speaking, GPSRs are less susceptible to the amount of excitation sources, but the proposed GPRLMs algorithms yielded better solutions than GPSRs in terms of spatial distribution under different excitation conditions, as shown in Fig. 7(a). From Fig. 7(c), the time cost raises with the increase of excitation sources, as the increase of measurement lead to the raised computational expense.

Fig. 7

Comparison results with different number of excitation nodes: (a) CNR, (b) LE, and (c) time.

JBO_22_4_045009_f007.png

To test the stability of the algorithm against noise, we added random Gaussian noise to the simulated data to mimic the noise caused by measurement error, registration error, and autofluorescence. For each simulation experiment, 10 independent runs were performed, because the amplitude of the simulated noise was randomly determined. Figure 8 shows the variation of LE and CNR (mean with standard deviation) with noise ranging from 5% to 25%. The simulation results indicate that the location accuracy of the target center provided by all the tested algorithms has no significant fluctuation (<0.03  mm) under different noise. The fluctuation range of CNR caused by noise is within 1 when the noise level ranges from 5% to 25%. In general, the tested algorithms perform quite stably in terms of LE and CNR when noise level is under 25%.

Fig. 8

(a) Variation of CNR (mean with standard deviation) with noise. (b) Variation of LE (mean with standard deviation) with noise.

JBO_22_4_045009_f008.png

Overall, the proposed GPRLM algorithms for joint 1 and Laplacian manifold regularization model performed quite stable under different excitation conditions and different noise levels.

3.3.

Reconstruction of Big Target

To test the big target reconstructed ability of the proposed algorithms, we conducted a simulation to reconstruct a big cylindrical target located at (12.8,11.5,16) mm, with a diameter of 4 mm and a height of 2 mm. The mouse model was identical to the setting in Sec. 3.2. As shown in Fig. 9, the results by GPRLMs have better consistency with the real target and are better localized around the corresponding target area while the reconstructions by GPSRs cannot fill the target region. The quantitative results listed in Table 3 also indicate the localization capability of GPRLMs is better than GPSRs.

Fig. 9

Reconstruction results of big target. (a), (b), (c), and (d) are cross-section views of the results obtained by GPSR-basic, GPRLM-basic, GPSR-BB, and GPRLM-BB, respectively. (e)–(h) are the 3-D views of the reconstructed target corresponding to (a)–(d).

JBO_22_4_045009_f009.png

Table 3

Comparison results in double-target simulations.

MethodReconstruction center (mm)CNRLE (mm)Time (s)
GPSR-basic(12.53, 10.94, 16.87)8.760.9255.21
GPRLM-basic(13.02, 11.13, 16.10)12.140.4089.61
GPSR-BB(12.53, 10.91, 16.87)8.770.9239.67
GPRLM-BB(12.76, 11.50, 16.42)14.870.8239.54
Note: The best results are in bold.

3.4.

Reconstruction of Multiple-Target with Different Sizes

Multiple-target resolving ability is equally important to reconstruction algorithms. Hence, another group of simulations was conducted to reconstruct two targets with different sizes and 1 mm separation. The small target is 2 mm in diameter and 2 mm in height. The big target is 4 mm in diameter and 2 mm in height. The double target is located in the liver with the centers at (14,11,16) mm and (14,7,16) mm, respectively. Obviously, the double-target setting is more difficult than the single-target case. The cross-section views of the results at the plane of Z=16  mm, and the 3-D views of the reconstructed target are shown in Fig. 10. It can be observed that the results by GPRLMs have less spreading and are better localized around the corresponding true center, while the reconstructions by GPSRs cannot distinguish the two targets. To further compare the spatial distribution of the solutions, Fig. 11 shows the profiles of the normalized fluorescent yield on XY plane at Z=16  mm along the line of X=14  mm. Two peaks can be easily recognized in the results of GPRLMs, whereas only one peak exists in the results of GPSRs. Based on Fig. 11, we computed the full width of half maximum (FWHM) values to assess the spatial aggregation of the solutions, which are listed in Table 4. More quantitative information of double-target results are presented in Table 5. The reconstruction results shown in this section demonstrate that the proposed GPRLMs for joint regularization model have stronger resolving capability and better location accuracy compared with GPSRs for 1-norm regularization model.

Fig. 10

Reconstruction results of double target with 1-mm boundary distance. (a), (b), (c), and (d) are cross-section views of the results obtained by GPSR-basic, GPRLM-basic, GPSR-BB, and GPRLM-BB, respectively. (e)–(h) are the 3-D views of the reconstructed target corresponding to (a)–(d).

JBO_22_4_045009_f010.png

Fig. 11

Profiles of normalized fluorescent yield on XY plane Z=16  mm along the line of X=14  mm, which corresponds to the black dashed lines in Figs. 10(a)10(d). (a) Comparison between GPSR-basic and GPRLM-basic. (b) Comparison between GPSR-BB and GPRLM-BB.

JBO_22_4_045009_f011.png

Table 4

Comparison results of FWHM in double-target simulations with different size.

MethodReal center (y-coordinate/mm)Reconstruction center (y-coordinate/mm)Real FWHM (mm)Reconstruction FWHM (mm)
GPSR-basic7N/A2N/A
1111.0044.03
GPRLM-basic76.3221.90
1111.2044.08
GPSR-BB7N/A2N/A
1110.9544.05
GPRLM-BB76.3221.86
1111.2244.06

Table 5

Comparison results in double-target simulations with different size.

MethodReconstruction center (mm)CNRLE (mm)Time (s)
GPSR-basic(15.04, 8.32, 15.22)6.321.85153.93
(13.31, 11.19, 16.58)0.92
GPRLM-basic(12.98, 7.46, 15.16)9.921.40172.77
(13.13, 11.42, 16.38)1.04
GPSR-BB(15.06, 8.33, 15.28)5.631.85137.25
(13.30, 11.15, 16.59)0.93
GPRLM-BB(13.01, 7.36, 15.10)8.651.39137.44
(13.30, 11.02, 16.62)0.94
Note: The best results are in bold.

3.5.

Validation with In Vivo Experimental Data

In this section, we further validate the proposed reconstruction method with in vivo experimental data. We used an integrated FMT/micro-CT dual-modality imaging system12 to acquire data. This system provides a stabilized compact red laser with wavelength at 670 mm and allows light transmission at the emission wavelength of 710 mm. Detailed configuration information of the imaging system can refer to Yi’s12 work. The animal studies were performed in accordance with the Fourth Military Medical University Guide for the Care and Use of Laboratory Animals formulated by the National Society for Medical Research.

In this experiment, a glass tube with 1.2-mm inner radius and 7.5-mm height was implanted into abdomen of an anesthetized mouse. The tube was injected with Cy5.5 solution (with the extinction coefficient of about 0.019  mm1μM1 and quantum efficiency of 0.23 at the peak excitation wavelength of 671 nm) to serve as the fluorescent target. According to the published paper,42 the fluorescent yield of Cy5.5 is 0.0402  mm1. From the structural data scanned by micro-CT, the actual center of the target was determined at (20.2,28.5,8.8) mm. Four excitation sources, we used are located at (20.16,32.32,10.88) mm, (20.16,11.04,10.88) mm, (26.88,21.76,10.88) mm, and (10.24,21.76,10.88) mm, respectively. For subsequent reconstruction, the CT data were discretized and segmented into five components, including heart, lungs, liver, kidneys, and muscle. Table 612,15,43 shows the main optical parameters about the excitation and emission processes used in reconstruction.

Table 6

Optical parameters of mouse organs at wavelength of 670 and 710 nm.

Organ670 nm710 nmg
μax(mm−1)μsx′(mm−1)μam(mm−1)μsm′(mm−1)
Muscle0.0750.4120.0430.3500.9
Heart0.0510.9440.0300.8700.85
Lungs0.1702.1570.0972.0930.94
Livers0.3040.6680.1760.6290.9
Kidneys0.0582.2040.0342.0210.86
Stomach0.0101.4170.0071.3400.92

Figure 12 shows the reconstruction results of the experimental data. The detailed evaluation results are listed in Table 7. Similar to the numerical studies presented in Sec. 3.2, GPRLMs performed better than GPSRs in terms of location accuracy and spatial aggregation. The LE by GPRLMs was around 0.8 mm, which was much smaller than that of GPSRs. As observed from Fig. 12, the results by our proposed GPRLMs for joint regularization model have fewer artifacts and are better localized around the glass tube.

Fig. 12

Reconstruction results of in vivo experiments. (a), (b), (c), and (d) are the fluorescent yield of reconstructed target at slice Z=8.4  mm with CT image, obtained by GPSR-basic, GPRLM-basic, GPSR-BB, and GPRLM-BB, respectively. The green circles in Figs. 9(a)9(d) indicate the location of real target. (e)–(h) are the fluorescent yield distribution of reconstructed target in 3-D views with CT image, obtained by GPSRs and GPRLMs, respectively.

JBO_22_4_045009_f012.png

Table 7

Evaluation results in in vivo experiments.

MethodReconstruction center (mm)CNRLE (mm)Time (s)
GPSR-basic(20.50, 29.79, 9.04)4.851.35177.61
GPRLM-basic(20.26, 29.29, 8.64)5.440.81149.11
GPSR-BB(20.60, 29.61, 8.94)7.281.1733.32
GPRLM-BB(20.48, 29.26, 8.67)8.230.8131.90
Note: The best results are in bold.

4.

Discussion and Conclusion

In this paper, we present an inverse model with joint 1 and Laplacian manifold regularization for FMT reconstruction. By incorporating prior information of spatial aggregation and sparsity of the fluorescent target observed in many FMT scenarios, the new joint model is expected to further improve the reconstruction quality. Compared with 1-norm regularization model, the proposed joint model has advantages in reconstruction of big size target or resolving multiple targets. Compared with TV penalty, the joint model is more geared to the features of distribution of tumor tissue and is easier to be applied in FMT modalities. Moreover, we also present two algorithms to solve the new inverse model based on gradient projection techniques. The effectiveness of the proposed reconstruction model and algorithms are validated with numerical studies and in vivo experimental data.

Simulation and experimental results demonstrate that better solution (fewer artifacts and smaller LE) is obtained by GPRLM due to the fine properties of manifold regularization. Numerical studies also indicate that the proposed reconstruction algorithms (GPRLM-basic and GPRLM-BB) performed quite stable with respect to measurement noise and change of excitation conditions. Compared with the widely used 1 regularization model, the proposed reconstruction method possesses a stronger ability of aggregation, which is well demonstrated.

The proposed joint model has several advantages over stand-alone TV or 1-norm regularization model. First of all, the complexity of Laplacian regularization is less than TV-norm. Second, Laplacian regularizer is easy to implement on FEM mesh. And last, the joint model gives consideration to both sparsity and spatial aggregation, which make it promising to apt to changing source distributions in different FMT applications.

The LE metric used in this paper is to assess the location accuracy of the center of mass, which cannot adequately reflect the reconstructed shapes. We recognize that target shapes are equally important as center localization. Although the joint model incorporating prior information of spatial aggregation and sparsity, the performance of the proposed algorithms (GPRLMs) was less than satisfactory in shape reconstruction. However, it should be pointed out that the proposed GPRLMs are preliminary methods for solving the joint 1 and Laplacian manifold regularization model. With a few modifications, many other sophisticated sparse recovery algorithms proposed in the field of compressed sensing can be used to solve the joint model and improve the reconstruction. Moreover, a more accurate light transport model such as the SP3 approximation may provide better performance in shape fitting. A future extension of this research can be the development of new efficient algorithms for solving this joint model and the corresponding performance assessment under different types of source distributions.

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NSFC) (61372046, 61401264, 11571012, and 61640418), the Science and Technology Plan Program in Shaanxi Province of China (2015KW-002), the Natural Science Research Plan Program in Shaanxi Province of China (2015JM6322), the Fundamental Research Funds for the Central Universities (GK201603025), Scientific Research Program Funded by Shaanxi Provincial Education Department under Grant No. (16JK1772), the Scientific Research Foundation of Northwest University under Grant Nos. (338050018 and 338020012) and the Project funded by China Postdoctoral Science Foundation under Grant No. (2016M602851). The authors would like to thank the School of Life Science and Technology of Xidian University for providing experimental data acquisition system.

References

1. 

R. Weissleder and M. Pittet, “Imaging in the era of molecular oncology,” Nature, 452 580 –589 (2008). http://dx.doi.org/10.1038/nature06917 Google Scholar

2. 

A. Ale et al., “FMT-XCT: in vivo animal studies with hybrid fluorescence molecular tomography-x-ray computed tomography,” Nat. Methods, 9 615 –620 (2012). http://dx.doi.org/10.1038/nmeth.2014 1548-7091 Google Scholar

3. 

D. Doleshel et al., “Targeted near-infrared imaging of the erythropoietin receptor in human lung cancer xenografts,” J. Nucl. Med., 53 304 –311 (2012). http://dx.doi.org/10.2967/jnumed.111.091124 JNMEAQ 0161-5505 Google Scholar

4. 

J. K. Willmann et al., “Molecular imaging in drug development,” Nat. Rev. Drug Discov., 7 591 –607 (2008). http://dx.doi.org/10.1038/nrd2290 NRDDAG 1474-1776 Google Scholar

5. 

D. T. Delpy and M. Cope, “Quantification in tissue near-infrared spectroscopy,” Philos. Trans. R. Soc. London Ser. B, 352 649 –659 (1997). http://dx.doi.org/10.1098/rstb.1997.0046 Google Scholar

6. 

V. Ntziachristos, J. Ripoll and R. Weissleder, “Would near-infrared fluorescence signals propagate through large human organs for clinical studies?,” Opt. Lett., 27 333 –335 (2002). http://dx.doi.org/10.1364/OL.27.000333 OPLEDP 0146-9592 Google Scholar

7. 

X. Zhang, C. Badea and G. Johnson, “Three-dimensional reconstruction in free-space whole-body fluorescence tomography of mice using optically reconstructed surface and atlas anatomy,” J. Biomed. Opt., 14 064010 (2009). http://dx.doi.org/10.1117/1.3258836 JBOPFO 1083-3668 Google Scholar

8. 

A. Tikhonov et al., Numerical Methods for the Solution of Ill-Posed Problem, Springer, Netherlands (2013). Google Scholar

9. 

D. Donoho et al., “Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit,” IEEE Trans Inf. Theory, 58 1094 –1121 (2012). http://dx.doi.org/10.1109/TIT.2011.2173241 Google Scholar

10. 

J. Ye et al., “Fast and robust reconstruction for fluorescence molecular tomography via a sparsity adaptive subspace pursuit method,” Biomed. Opt. Express, 5 387 –406 (2014). http://dx.doi.org/10.1364/BOE.5.000387 BOEICL 2156-7085 Google Scholar

11. 

J. Baritaux, K. Hassler and M. Unser, “An efficient numerical method for general lp regularization in fluorescence molecular tomography,” IEEE Trans. Med. Imaging, 29 1075 –1087 (2010). http://dx.doi.org/10.1109/TMI.2010.2042814 ITMID4 0278-0062 Google Scholar

12. 

H. Yi et al., “Reconstruction algorithms based on l1-norm and l2-norm for two imaging models of fluorescence molecular tomography: a comparative study,” J. Biomed. Opt., 18 056013 (2013). http://dx.doi.org/10.1117/1.JBO.18.5.056013 JBOPFO 1083-3668 Google Scholar

13. 

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

14. 

S. J. Kim et al., “An interior point method for large-scale l1-refularized least squares,” IEEE J. Selected Topics in Signal Process., 4 606 –617 (2007). http://dx.doi.org/10.1109/JSTSP.2007.910971 Google Scholar

15. 

H. Guo et al., “Improved sparse reconstruction for fluorescence molecular tomography with l1/2 regularization,” Biomed. Opt. Express, 6 1648 –1664 (2015). http://dx.doi.org/10.1364/BOE.6.001648 BOEICL 2156-7085 Google Scholar

16. 

D. Zhu et al., “Comparison of regularization methods in fluorescence molecular tomography,” Photonics, 1 95 –109 (2014). http://dx.doi.org/10.3390/photonics1020095 Google Scholar

17. 

D. Han et al., “Efficient reconstruction method for l1 regularization in fluorescence molecular tomography,” Appl. Opt., 49 6930 –6937 (2010). http://dx.doi.org/10.1364/AO.49.006930 APOPAI 0003-6935 Google Scholar

18. 

Z. Xu et al., “l1/2 regularization,” Sci. China Inf. Sci., 53 1159 –1169 (2010). http://dx.doi.org/10.1007/s11432-010-0090-0 Google Scholar

19. 

J. Dutta et al., “Joint l1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol., 57 1459 –1476 (2012). http://dx.doi.org/10.1088/0031-9155/57/6/1459 PHMBA7 0031-9155 Google Scholar

20. 

J. Feng et al., “Total variation regularization for bioluminescence tomography with the split Bregman method,” Appl. Opt., 51 4501 –4512 (2012). http://dx.doi.org/10.1364/AO.51.004501 APOPAI 0003-6935 Google Scholar

21. 

M. Freiberger, C. Clason and H. Scharfetter, “Total variation regularization for nonlinear fluorescence tomography with an augmented Lagrangian splitting approach,” Appl. Opt., 49 3741 –3747 (2010). http://dx.doi.org/10.1364/AO.49.003741 APOPAI 0003-6935 Google Scholar

22. 

C. Li et al., “An efficient augmented Lagrangian method with applications to total variation minimization,” Comput. Optim. Appl., 56 507 –530 (2013). http://dx.doi.org/10.1007/s10589-013-9576-1 CPPPEF 0926-6003 Google Scholar

23. 

S. C. Davis et al., “Image-guided diffuse optical fluorescence tomography,” Opt. Exp., 15 4066 –4082 (2007). http://dx.doi.org/10.1364/OE.15.004066 OPEXFF 1094-4087 Google Scholar

24. 

W. He et al., “Surface fluorescence molecular tomography with prior information,” Appl. Opt., 53 402 –409 (2014). http://dx.doi.org/10.1364/AO.53.000402 APOPAI 0003-6935 Google Scholar

25. 

R. Holt, S. Davis and B. Pogue, “Regularization functional semi-automated incorporation of anatomical prior information in image-guided fluorescence tomography,” Opt. Lett., 38 2407 (2013). http://dx.doi.org/10.1364/OL.38.002407 OPLEDP 0146-9592 Google Scholar

26. 

A. L. Vahrmeijer et al., “Image-guided cancer surgery using near-infrared fluorescence,” Nat. Rev. Clin. Oncol., 10 507 –518 (2013). http://dx.doi.org/10.1038/nrclinonc.2013.123 Google Scholar

27. 

J. Tenenbaum, V. Silva and J. Langford, “A global geometric framework for nonlinear dimensionality reduction,” Science, 290 2319 –2323 (2000). http://dx.doi.org/10.1126/science.290.5500.2319 SCIEAS 0036-8075 Google Scholar

28. 

S. Roweis and L. Saul, “Nonlinear dimensionality reduction by locally linear embedding,” Science, 290 2323 –2326 (2000). http://dx.doi.org/10.1126/science.290.5500.2323 SCIEAS 0036-8075 Google Scholar

29. 

Z. Wang, N. Nasrabadi and T. Huang, “Semisupervised hyperspectral classification using task-driven dictionary learning with Laplacian regularization,” IEEE Trans. Geosci. Remote Sens., 53 1161 –1173 (2015). http://dx.doi.org/10.1109/TGRS.2014.2335177 IGRSD2 0196-2892 Google Scholar

30. 

A. Babaeian et al., “Angle constrained path for clustering of multiple manifolds,” in IEEE Int. Conf. Image Processing, 4446 –4450 (2015). http://dx.doi.org/10.1109/ICIP.2015.7351647 Google Scholar

31. 

W. Bian et al., “Minimizing nearest neighbor classification error for nonparametric dimension reduction,” IEEE Trans. Neural Network Learn. Syst., 25 1588 –1594 (2014). http://dx.doi.org/10.1109/TNNLS.2013.2294547 Google Scholar

32. 

M. Belkin and P. Niyogi, “Laplacian eigenmaps for dimensionality reduction and data representation,” Neural Comput., 15 1373 –1396 (2013). http://dx.doi.org/10.1162/089976603321780317 NEUCEB 0899-7667 Google Scholar

33. 

M. Figueiredo, R. Nowak and S. Wright, “Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems,” IEEE J. Sel. Top. Signal Process., 1 586 –597 (2007). http://dx.doi.org/10.1109/JSTSP.2007.910281 Google Scholar

34. 

A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express, 13 9847 –9857 (2005). http://dx.doi.org/10.1364/OPEX.13.009847 OPEXFF 1094-4087 Google Scholar

35. 

V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng., 8 1 –33 (2006). http://dx.doi.org/10.1146/annurev.bioeng.8.061505.095831 ARBEF7 1523-9829 Google Scholar

36. 

X. Song et al., “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express, 15 18300 –18317 (2007). http://dx.doi.org/10.1364/OE.15.018300 OPEXFF 1094-4087 Google Scholar

37. 

D. Wang, Y. C. X. Liu and J. Bai, “A novel finite-element-based algorithm for fluorescence molecular tomography of heterogeneous media,” IEEE Trans. Inf. Technol. Biomed., 13 766 –773 (2009). http://dx.doi.org/10.1109/TITB.2009.2015144 Google Scholar

38. 

L. Zhao et al., “lp regularization for early gate fluorescence molecular tomography,” Opt. Lett., 39 4156 –4159 (2014). http://dx.doi.org/10.1364/OL.39.004156 OPLEDP 0146-9592 Google Scholar

39. 

X. Song et al., “Automated region detection based on the contrast-to-noise ratio in near-infrared tomography,” Appl. Opt., 43 1053 –1062 (2004). http://dx.doi.org/10.1364/AO.43.001053 APOPAI 0003-6935 Google Scholar

40. 

B. Dogdas et al., “Digimouse: a 3D whole body mouse atlas from ct and cryosection data,” Phys. Med. Biol., 52 577 –587 (2007). http://dx.doi.org/10.1088/0031-9155/52/3/003 PHMBA7 0031-9155 Google Scholar

41. 

H. Guo et al., “Adaptive hp finite element method for fluorescence molecular tomography with simplified spherical harmonics approximation,” J. Innov. Opt. Health Sci., 7 1350057 (2014). http://dx.doi.org/10.1142/S1793545813500570 Google Scholar

42. 

M. A. Naser and M. S. Patterson, “Improved bioluminescence and fluorescence reconstruction algorithms using diffuse optical tomography, normalized data, and optimized selection of the permissible source region,” Biomed. Opt. Express, 2 169 –184 (2011). http://dx.doi.org/10.1364/BOE.2.000169 BOEICL 2156-7085 Google Scholar

43. 

G. Alexandrakis, F. 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

Biography

Xuelei He received his BS degree from Wuhan University in 2015. He is a PhD student in the School of Information Sciences and Technology, Northwest University, Xi’an, China. His current research interests include tomography and optimization.

Xiaodong Wang received his BS degree from Harbin Institute of Technology, Harbin, China, in 1998 and his MS degree from Inner Mongolia University of Technology, Hohhot, China, in 2007, and his PhD in computer application technology from Xidian University, Xi’an, China, in 2014. He is a postdoctoral research at the Information Science and Technology Institute, Northwest University, Xi’an, China. His current research interests include convex optimization, tomography, and pattern recognition.

Huangjian Yi received her BS degree in applied mathematics and her PhD in intelligent information processing from Xidian University, Xi’an, China, in 2008 and 2013, respectively. She is currently a lecturer in the School of Information Sciences and Technology, Northwest University, Xi’an, China. Her research interests include optics, molecular imaging, medical image processing, and the applications of molecular imaging.

Yanrong Chen received her BS degree from Central South University in 2014. She is a PhD student in the School of Information Sciences and Technology, Northwest University, Xi’an, China. Her current research interest is image segmentation.

Xu Zhang received his BS degree from Northwest University in 2014. He is a master’s student in the School of Information Sciences and Technology, Northwest University. His current research interests include tomography and image processing.

Jingjing Yu received her MS degree from Xi’an Jiaotong University in 2005 and her PhD from Xidian University in 2012. She is an assistant professor in the School of Physics and Information Technology, Shaanxi Normal University, Xi’an, China. Her current research interest is intelligent information processing.

Xiaowei He received his MS degree from Xi’an Jiaotong University in 2005 and his PhD from Xidian University in 2011. He is a professor in the School of Information Sciences and Technology, Northwest University, Xi’an, China. His current research interest is optical molecular imaging.

© 2017 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2017/$25.00 © 2017 SPIE
Xuelei He, Xiaodong Wang, Huangjian Yi, Yanrong Chen, Xu Zhang, Jingjing Yu, and Xiaowei He "Laplacian manifold regularization method for fluorescence molecular tomography," Journal of Biomedical Optics 22(4), 045009 (21 April 2017). https://doi.org/10.1117/1.JBO.22.4.045009
Received: 7 December 2016; Accepted: 7 April 2017; Published: 21 April 2017
Lens.org Logo
CITATIONS
Cited by 18 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Reconstruction algorithms

Luminescence

Tomography

Fluorescence tomography

3D modeling

Mouse models

Computer simulations

Back to Top