bioluminescence tomography ( blt ) is a rapidly growing field of research in optical molecular imaging , which allows for the visualization of normal and abnormal cellular processes in living subjects at the molecular or genetic level [ 14 ] .
with blt , we seek to recover the spatial distribution of bioluminescent light source inside a small animal from external noninvasive measurements . generally speaking
, the internal source intensity is closely related to the strength of the molecular / cellular activity , such as gene expression .
thus , this imaging modality can provide in - depth information of the internal biological sources concerned in longitudinal monitoring and quantitative assessment changes and efficacy and thus further facilitates our understanding of bio - molecular processes as they occur in living animals . when using blt technique to measure efficiency of a genic therapy or to observe the growth or migration of cancer cells , accurate detection of different sources that differ greatly in density or power is instrumental
; for example , it may yield a great deal of information regarding tumor dissemination and burden in various sites before the development of gross disease [ 1 , 7 , 8 ] .
therefore , the emphasis of this paper is multiple - source reconstruction that has not been sufficiently considered to date in blt .
most reconstruction methods for blt can be classified to model - based reconstruction . in this case ,
given a light propagation model , the flux on the boundary can be predicted with numerical methods such as the finite element method ( fem ) by combing with the structural information and optical parameters regarding different organs .
and then the blt is formulated as an optimization problem of minimizing the discrepancy between the boundary measurements and the predicted light intensities on the tissue surface . in the reconstruction procedure , the ill posedness of the blt problem does pose a challenge for determining a unique solution of the tomographic problem .
different strategies have been proposed for coping with the ill posedness of blt inverse problems .
these studies obtain stable reconstruction by increasing the amount of independent measurements with spectrally resolved approaches [ 1113 ] , or by reducing the number of unknowns [ 10 , 14 ] , or with regularization techniques to incorporate some a   priori information regarding the inverse source problem [ 1517 ] . in this paper
, we focus our attention on the multiple - source reconstruction with monochromatic boundary measurements where regularization techniques are inevitable in the reconstruction process .
the existing regularization - based reconstruction schemes in bioluminescent imaging to date can be loosely classified into three categories : l2 regularization , l1 regularization , and implicit regularization such as tsvd and lsqr [ 18 , 19 ] . through regularization ,
some constraints are applied to reconstruction and yield an approximate solution of the blt problem .
no matter which regularizer is used , source location and visualization are still needed for preclinical practice .
most source location schemes are directly based on the reconstructed density vector and the larger the density , the more probable the source center . specifically , according to a priori knowledge of the number of sources , several nodes with larger density values
are identified as the promising sources or set a global threshold by referring to the maximum density and only those nodes with a density value higher than the threshold will be displayed . in most applications of blt , for example ,
monitoring cancer metastasis , neither the sources number nor an appropriate global threshold is easy to determine .
this is mainly due to the fact that bioluminescent lights are usually weak and diffuse , and consequently the number of potential sources is hard to estimate only by surface photon distributions .
moreover , the global threshold strategy is unfeasible for distinguishing multiple sources with distinct difference in power . especially in l2 norm regularization cases ,
the obtained solution is usually oversmoothing , and thus a lower threshold will incur some artifacts in the final images whereas a higher one will discard some small potential sources . consequently , effective reconstruction scheme for multiple sources with different powers deserves further investigation . in this paper , we develop a finite element mesh aggregating approach for multiple - source reconstruction in blt .
first , we propose a multiple - source detecting strategy . rather than assuming independence between mesh nodes , the proposed reconstruction strategy exploits spatial structure of the nodes and characteristic of energy decay to adaptively
second , we integrate the proposed reconstruction strategy with regularization - based inverse algorithms to build a unified framework for solving blt inverse problem . numerical simulations and phantom experiments demonstrate the effectiveness of this framework . the paper is organized as follows . in section 2
, we present a multiple - source reconstruction framework with the emphasis on the finite - element - mesh - aggregating - based source detection strategy . in section 3
section 4 presents a phantom experiment to further test the effectiveness of the proposed method .
radiative transfer equation ( rte ) plays an important role in image reconstruction by predicting the bioluminescence light intensities on the tissue boundary , but solving rte remains an intractable task for biological tissue with spatially nonuniform optical properties and complex tissue geometries . instead , some approximations to rte have been established to overcome the difficulty of directly solving rte . among them , the diffusion approximation ( da ) model has been extensively used to describe the photon propagation in tissue where there is scattering dominant absorption [ 514 ] . here ,
we restrict our discussion to the da model for simplicity . the steady state diffusion equation complemented with the robin boundary condition
can be expressed as follows : 
 ( 1)(d(r)(r))+a(r)(r)=s(r ) ,  ( r ) , 
 ( 2)(r)+2a(r;n , n)d(r)(v(r)(r))=0 ,  ( r ) , 
							 where (r ) is the photon power density at r   , s(r ) is an isotropic source distribution of gene expression , and d(r ) and a(r ) are the optical diffusion and absorption coefficient , respectively . in this work
the term v(r ) in ( 2 ) denotes the unit outer normal at boundary  , a(r ; n , n )  ( 1 + r(r))/(1  r(r ) ) is the boundary mismatch factor accounting for different refractive indices across the boundary . following the standard finite element analysis , support domain  is discretized into t vertex nodes ( n1 , n2 ,  , nt ) and ne mesh elements , denoted as   ( l = 1,2 ,  , ne ) ; then (r ) and source term s(r ) can be approximately expressed as 
 ( 3)(r)h(r)=k=1tkk(r ) ,  r,s(r)sh(r)=k=1tskk(r ) ,  r , 
							 where k is the approximate nodal value of (r ) on the kth node nk , k(r ) the nodal basis function with support over the elements  , sk the discretized nodal values of s(r ) , and k(r ) the interpolation basis functions , which is usually the same with k(r ) . based on ( 1)(3 ) , a matrix equation of the linear relationship between source distribution and boundary measurements can be derived [ 10 , section 2 ] : 
 ( 4)as= , 
							 where a is a typical ill - conditioned matrix and  * represents measurable boundary nodal photon density .
in real blt experiments ,  * is computed from the surface flux image captured with a ccd camera . as mentioned in section 1
, the flux density on the boundary can be predicted according to a forward model , thereby a natural choice for source reconstruction is to minimize the misfit between predicted data and measurements , that is , 
 ( 5)s = arg  mins||as||2 . 
							 to deal with the ill posedness of blt inverse problem , permissible source region
is usually incorporated into the reconstruction model by spatially constraining the reconstruction domain to the area of interest [ 10 , 14 , 16 , 23 ] .
a more effective approach to reconstruction is using regularization to act as an algebraic stabilizer in estimating solutions . using a general lp  ( 0 < p  2 ) norm constraint , we reformulate the objective function for blt reconstruction in ( 5 ) : 
								 ( 6)sreg = arg  mins{||as||22+||s||p } , 
							 where the first term represents reconstruction error and the second is regularization term that fuses a priori knowledge or constrains into reconstruction .
regularization parameter  > 0 provides a tradeoff between data fitting and constraints regarding solutions .
obviously , tikhonov regularization method is a special case of ( 6 ) for p = 2 , that is , using an l2-norm regularizer . for p = 1 ,
l1-norm - based sparse regularization methods have recently attracted considerable amount of attention in blt [ 17 , 2325 ] and the reconstructions results therein witnessed some improvements in image quality .
based on the solution ( a source density vector ) obtained in section 2.2 , source localization and imaging is then performed by combining with fem mesh information . facing the dilemma of threshold choice mentioned in section 1
, we are hoping for an adaptive method that can avoid the difficulty of threshold selection while at the same time removing artifacts in the reconstructed images with relatively lower computational cost .  
consider that in most applications of blt , for example , detecting events that occur during the early stages of disease progression , the bioluminescent sources we want to recover are often localized in some small subregions of the domain .
on the other hand , because light intensity is heavily attenuated in biological tissue and falls off exponentially from the illumination point , the diffusion range of a bioluminescent source is limited by the source strength .
consequently , when taking the spatial structure of the mesh nodes into account , the source density vector should have a spatial aggregation on the mesh , which is also illustrated in the experiments in section 3 ( figure 4 ) .
it is found that , in a very small local region , if a node in the mesh has a maximum density value , with a very high probability its adjacent nodes are also with larger density .
it is found that in a very small local region , if a node in the mesh has a maximum density value , with a very high probability its adjacent nodes are also with a larger density .
we also observe that there are some nodes with smaller density in the vicinity of nodes with the larger density .
these observations are helpful for discriminating pseudosource from a cluster of mesh nodes and removing artifacts in images . on the basis of the above analysis , an iterative multiple - source detection strategy ( msds )
consequently , the preprocessing of solution with a small threshold of cmax  ( si ) is helpful to remove pseudosources and reduce the data size to be processed in the subsequent steps .
consequently , the preprocessing of solution with a small threshold of cmax  ( si ) is helpful to remove pseudosources and reduce the data size to be processed in the subsequent steps .
= { si | i  n , si > 0}. 
 step 4initial the sources number k = 1
we can find out the other elements that directly adjoin the node j , if any , according to the mesh structure information .
we can find out the other elements that directly adjoin the node j , if any , according to the mesh structure information .
step 6if set o is null , stop ; otherwise k :   = k + 1 , and go to step 5 . 
 
if set o is null , stop ; otherwise k :   = k + 1 , and go to step 5 . with the steps defined above , we provide an automatic method to estimate the number of sources from the reconstruction results iteratively .
the final results contain k sources . here , k is the number of subsets of the initial set o obtained at the end of the above iteration .
when pi  ( i = 1 , l , k ) has more than one member , we call this situation  overrepresentation ,  the nodes related to these elements will aggregate to represent a single source and the node with largest density value sj is regarded as the source center for simplicity .
eventually , the cartesian coordinates of the reconstructed sources are obtained by their node index in the finite element mesh .
based on the foregoing reconstruction scheme , we build a unified regularization framework for multiple - source reconstruction by integrating the msds with the general lp - norm regularization , as shown in figure 1 .
the msds is a relatively independent component of the framework , and hence different regularizer and different reconstruction algorithms can be utilized according to the practice of blt .
in this section , we present some numerical experiments to demonstrate the utility and the effectiveness of the proposed method in multiple - source settings . comparison is performed between the proposed msds and the traditional global threshold strategy ( gts ) .
it should be pointed that the main theme of this paper is to evaluate the performance of this framework for multiple - source reconstruction in blt , rather than the comparison between specific reconstruction algorithms . as representatives of algorithms using l1 and l2 regularization , tikhonov regularization method   and l1ls   are , respectively , combined with the above two strategies to recover the interior source distribution from the synthetically boundary measurements .
consequently , the reconstruction methods evaluated in the following experiments include tikhonov + msds , tikhonov + gts , l1ls + msds , and l1ls + gts .
it is known that regularization parameter is crucial to yield a good solution for ill - posed problems , and the choice of regularization parameter is usually nontrivial . in this paper , the regularization parameter for tikhonov method was determined with the adaptive method proposed in . as for l1ls , the parameter  was chosen as suggested in , that is ,  = 0.1||2a*||. all the experiments were performed on a cylindrical mouse chest numerical phantom as shown in figure 2(a ) .
the specific optical properties of different organs are listed in table 1 . in the first study
two sphere sources with radius of 0.5  mm were positioned in the left lung with the centers at s1 = ( 9 , 3.5,15 ) and s2 = ( 9,3.5,15 ) , respectively .
they were uniform in size and shape . to illustrate the point of our discussion ,
we consider four cases of experiment settings : ( i ) both of the initial source densities were 1  nw / mm ; ( ii ) to ( iv ) the densities of s1 were still 1  nw / mm , but the densities of s2 were 0.5  nw / mm , 0.25  nw / mm , and 0.125  nw / mm , respectively , that is , the ratios of the power of source s2 to that of s1 were 2  :  1 , 4  :  1 , and 8  :  1 . in the following experiments ,
the model was discretized into a fine tetrahedral element mesh and synthetic measurements were generated by solving the forward model with fem . to simulate the noise involved in real blt experiment ,
figures 2(b)2(e ) show the forward mesh and the simulated photon distribution on the surface in the above four source settings . obviously , it is difficult to predict the source number only according to the photon distribution especially in case ( iii ) and case ( iv ) . in the reconstruction process ,
a permissible source region strategy was also employed as a priori information to decrease the ill posedness of blt inverse problem , which was defined as { ( x , y , z ) | 8 < ( x + y ) < 12,13.5 < z < 16.5 } .
following the proposed reconstruction framework the reconstructions were carried out with the aforementioned four methods under different source settings .
the first row and the third row of figure 3 show the final reconstruction results by tikhonov method and l1ls method combined with the proposed msds . for comparison , the second row and the fourth row of figure 3 present the corresponding reconstructed results rendered from gts , where a global threshold ( 35% of the maximum density value ) was used .
it is obvious that the two sources are accurately detected by the proposed msds combined with different regularization methods in all the cases considered . on the other hand , for case ( iii ) and case ( iv ) , only the source with larger power
is detected by tikhonov + gts and l1ls + gts , whereas the other weaker one is lost in the final reconstruction results .
to quantitativly assess reconstruction results in different power settings , we summarize location errors and reconstructed powers by different reconstruction schemes in table 2 , where the second column represents the actual initial power ratio of s1 to s2 , and s1 and s2 denote the corresponding reconstructed sources .
, it is seen that l1-norm - based method l1ls generally performs better than l2-norm - based tikhonov method in terms of reconstructed powers and locations . 

figure 4 illustrates the mesh aggregating process of msds and compares the final reconstruction results of msds with those of gts in case ( i ) .
we can observe that there are some nodes with smaller density value in the vicinity of the two nodes with larger density , as shown in figures  4(a ) and 4(b ) .
apparently , retaining all of the nonzero components of the regularized solution will incur some artifacts in the final reconstruction image , in particular for l2 norm solution by tikhonov regularization method .
the results in figures  4(c ) and 4(d ) show that the traditional gts directly discards those nodes with density value lower than the given threshold in the final results to improve the image quality .
usually , a higher threshold is preferred in the literature , thus a threshold of 0.35max  ( si ) was used in the experiments for gts method [ 16 , 29 ] . as a result , those suspect targets with density lower than threshold will be omitted in this way .
unlike traditional methods , the proposed msds considers not only density value of a node but also mesh structure used in reconstruction and thus it has an ability to remove pseudosources and retain weak suspect sources in the final reconstruction results , as shown in figures  4(e)-4(f ) and 3 . in the second experiment
, we attempt to reconstruct sources with synthetic data generated from four scattered sources with different initial powers , which may be a common case in tumor metastasis .
specifically , the power setup was according to ratio of 8  :  4  :  2  :  1 and the maximum power density was 1  nw / mm .  
figure 5 shows 3d views of the results of tikhonov regularization method and l1ls method , respectively , combined with gts and msds .
obviously , it is hard for traditional gts method to detect multiple sources with lower power density in such experimental setting , whereas the proposed msds accurately distinguishes all of the sources . in view of the idea that the proposed multiple - source reconstruction approach utilizes underlying mesh structure information ,
it is necessary to assess the influence of different fem discretization on the proposed method .
therefore , we conducted a set of double - source experiments under different discretization level .
the results in figure 6 ( where the number of nodes in reconstruction domain denotes different discretization level or mesh size ) show the influence of finite element mesh on reconstruction . for tikhonov regularization method combined with msds ,
the location error increases slightly after a decrease along with the increasing of mesh size and the reconstructed power presents a similar variation trend . as for l1ls combined with msds ,
finite element discretization does affect reconstructed results in the sense that the location error and the reconstructed power vary with the change of mesh . however , for all of the discretization levels considered , the proposed method is able to accurately localize and quantify light source distribution .
this set of blt experiments were conducted with a dual - modality blt / micro - ct system [ 17 , 30 ] .
a backthinned , backilluminated cooled ccd camera is used to measure the signal on the phantom surface from four directions at 90-degree intervals .
the heterogeneous mouse chest phantom with 30  mm height and 15  mm diameter consists of four parts that represent muscle , lungs , heart , and bone , respectively .
two small holes of diameter 2  mm were drilled in the phantom to place glass capillary with 1  mm inside diameter .
luminescent solutions of height 2  mm were extracted from a red luminescent light stick ( glow products , canada ) and then injected to glass capillary to serve as one testing source .
the real center positions of the two testing sources were ( 9,2 , 16.6 ) and ( 9 , 3,16.6 ) .
we collected 100 gray level images of the sources , which were taken by the ccd camera every one minute .  
figure 7 shows the fitted decay curve of light density . according to the decay curve
, we can obtain sources with different intensities by controlling the injection time of luminescent solutions .
three groups of experiments were conducted , and the ratios of the intensity of source s2 to that of s1 were 1  :  1 , 2  :  1 , and 4  :  1 , respectively . figures 8(a)8(c ) show the front views of the corresponding measured data on ccd under different intensity settings .
subsequently , a permissible source region was roughly determined according to the surface flux density distribution , which is expressed as { ( x , y , z ) | 8 < ( x + y ) < 13,15 < z
< 18}. the phantom model was discretized into 4202 nodes and 21721 tetrahedra . after mapping the collected optical data on the three - dimensional phantom surface
, we performed four rounds of reconstruction with tikhonov + gts , l1ls + gts , tikhonov + msds , and l1ls + msds under different source intensity settings .
the normalized reconstruction results of tikhonov regularization method are similar to that of l1ls . to avoid interminable description , figure 9 only presents comparison results between tikhonov + gts and tikhonov + msds . for all of the testing cases considered in phantom experiments , tikhonov + msds and l1ls + msds can accurately detect two sources , and the maximum location error is 1.7  mm . even for the case of real intensity ratio 4  :  1 ,
the reconstructed source strength ratios of them were 3.12  :  1 and 2.97  :  1 .
in stark contrast to the proposed methods , traditional global threshold methods failed to reconstruct the weaker of the two sources , as shown in figure 9(c ) . compared with the results of using gts ( the top row of figure 9 ) , the proposed msds methods produce fewer artifacts in the reconstructed images ( the bottom row of figure 9 )
accurately reconstructing and distinguishing several sources with different intensities is a challenge problem in blt , which is also an essential ability for serial observation of disease progression or response to therapy in the same animal over time . in this work
, we present a unified framework for multiple - source reconstruction by integrating a novel multiple - source detection strategy with regularization - based reconstruction process .
the effectiveness of this regularization framework is validated with numerical simulations and further confirmed with phantom experiments .
first , there is no need for a prior knowledge regarding source number , which is automatically estimated from the reconstruction results iteratively .
second , the regularization framework is general since it can work with different regularizers and inverse algorithms .
the proposed msds is also easily applied to other finite - element - based reconstruction schemes to improve the final reconstruction results or image quality .
there are several limitations to the proposed method . as indicated in the experiment results , sparseness - inducing regularization method ( l1ls ) performs better than l2 norm method ( tikhonov ) .
this is mainly because l1 norm solution accords with the sparsity nature of bioluminescent source distribution in these applications .
consequently , how to select appropriate regularizer and inverse algorithm for specific blt application is very important when using this framework .
lp(0 < p < 1 ) norm regularized reconstruction has been tried for recovery of signals with weak sparsity in other image processing fields . so far , related researches have not yet been reported in blt . based on the proposed regularization framework
, our future studies will investigate the effectiveness of other forms of regularizer for the ill - posed inverse problem of blt .
although only the da model is considered for the sake of simplicity , the proposed blt reconstruction framework has no limitation on the forward model .
the performance of our framework might be improved by using more accurate forward models , which is also the direction of our further work .
 a finite element mesh aggregating approach is presented to reconstruct images of multiple internal bioluminescence sources . rather than assuming independence between mesh nodes , the proposed reconstruction strategy exploits spatial structure of nodes and aggregation feature of density distribution on the finite element mesh to adaptively 
 determine the number of sources and to improve the quality of reconstructed images . 
 with the proposed strategy integrated in the regularization - based reconstruction process , reconstruction algorithms need no a priori knowledge of source number ; even more importantly , they can automatically reconstruct multiple sources that differ greatly in density or power . 
1. Introduction
2. Multiple-Source Reconstruction Framework
3. Numerical Results and Analysis
4. Phantom Study
5. Discussions and Conclusion