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 .