
				  the cell morphodynamics underlying the processes that control the organized spatial distribution of cells during the embryonic development are nowadays largely unknown .
a single fertilized egg - cell ( zygote ) undergoes cellular divisions , differentiations , and interactions giving rise to specific forms of tissues and organs .
the reconstruction of the nuclei shape and the identification of their number and position during embryogenesis are an essential tool for an integrated understanding of such morphogenetic processes .
this kind of information is relevant to estimate the cell growth and identify cell divisions and apoptosis .
it is a helpful tool to extract parameters such as the cell proliferation rate in time and space .
it is relevant for investigating stem cells populations and early steps of cancerogenesis , opening the way for the preclinical evaluation of anticancer drug effects in vivo .
we apply in this paper a region - based segmentation algorithm to segment the nuclei in a live zebrafish embryo .
since during the acquisition the embryo is growing up , the cell number strongly increases .
furthermore , since prior and during mitosis the nuclei condense and lengthen along the future cell division plane , the segmentation is a helpful tool to identify cell divisions .
we consider and compare two different forms of the equation originally proposed by chan and vese .
comparison is performed through qualitative and quantitative analysis of results on both phantoms and real data .
we then apply to our dataset the formulation that ensures the best performances . in the rest of the paper
we introduce the active contour model without edges and both tested formulations . in section 4
the zebrafish ( danio rerio ) is a model organism extensively studied to explore the gene function and vertebrate development [ 27 ] and might soon become a major model organism for preclinical drug testing by pharmaceutical industries .
although distant from humans , nevertheless it has comparable organs and tissues and may supplement higher vertebrate models . compared with mice , zebrafish embryo develops rapidly and externally to the mother , favoring its manipulation .
its transparency allows the visual inspection of the internal anatomy . due to its ability to regenerate fins , skin , the heart , and the brain
, it is an object of study to understand healing / repair mechanisms in vertebrates . 
				  to reconstruct the shape and
the position of every nucleus in a live zebrafish embryo is necessary to stain all the cells and use an acquisition technique with micrometrical resolution .
this is a very challenging approach requiring specific methodologies and tools in embryos engineering and microscopy imaging .
recent advances in imaging strategies open the way to in 4d imaging of live animals with a resolution at the cellular level and enough contrast to allow segmentation of individual cells . through the ubiquitous expression of fluorescent proteins in the zebrafish embryo ,
it is possible to label all the cells and perform time - lapse clsm ( confocal laser scanning microscopy ) imaging throughout embryonic development .
the embryo used in this paper to design and validate our strategy has been stained through injection at the one cell stage of rnas encoding farnesylated mcherry fluorescent protein and histone h2b / egfp fusion protein .
this allows the simultaneous acquisition of two different signals , membranes and nuclei , respectively . in this paper
the images have been acquired using a confocal microscope leica sp2 aobs with a 40x/0.8 na water objective and 488 nm and 561 nm laser light excitation .
we used the best compromise in terms of embryo survival and spatial and temporal resolution , as suggested in .
the data has been acquired for four hours ( 25c under the microscope ) , starting at 3.5-hour postfertilization ( development at 28c ) from the animal pole . due to the lower temperature , during the acquisition the normal embryo development  
is slowed down and by the end it reaches the shield stage ( see figure 1 ) .
images are acquired parallel to the focal plane xy , at regularly spaced z depth to build a z stack .
the procedure has been repeated every 5 minutes to generate a temporal sequence of 49 volumes .
the dimension of each voxel is 0.58  0.58  1.04 . in the entire period of development , the whole embryo looks like a sphere with a diameter of 740 . as the 3d images have a physical dimension of 300  300  31  , they cover only the top part of the embryo , as depicted in figure 1 . observing figure 2(a )
is possible to note how the typology of staining influences the features of original data . as the histone h2b protein is a component of the chromatin , nuclei in mitosis
conversely cell nuclei in interphase show a nonuniform image intensity and boundaries often not clearly defined .
we explain in next section how this peculiarity leads to the choice of the segmentation method . 
				

				  one of the most popular approaches for image segmentation is the active contour model or snake , based on the evolution of a curve attracted by image boundaries in order to detect objects .
implemented such models using level set methods , introduced first by osher and sethian   and extensively used to track the evolution of fronts in a variety of applications . with the level set methods the desired curve of active contour models
is embedded as zero level set of an implicit function , overcoming the drawbacks of the original snakes models   as the difficulties related to topological transformations .
however , the above - mentioned active contour models are edge - based ; that is , a curve is evolved with a speed function depending on a precomputed edge indicator defined by the image gradient .
dealing with image derivatives , a preprocessing step is required to remove the noise and smooth the image .
if the objects are not well contrasted , the denoising can disrupt image information and preclude the correct identification of objects surface . lately
a region - based level - set method has been proposed by chan and vese   to segment objects whose boundaries are not defined by a gradient .
the model does not require preprocessing and is based on the informations provided directly by the image intensity .
the chan vese approach is then well adapted to situations in which images are noisy , and it is difficult to extract the boundary of the desired object .
let us consider a 3d image i as a real positive function defined in a rectangular domain   . the chan vese   model moves an arbitrary surface s   for decomposing the image i into two regions of approximately piecewise - constant intensities : the objects to be detected and the background .
such decomposition is obtained by minimizing the following energy functional : 
								 ( 1)e(s , si , so)=(area(s ) )  + iinside(s)|i(x , y , z)si|2dx  dy  dz  + ooutside(s)|i(x , y , z)so|2dx  dy  dz , 
							 where si and so are , respectively , the averages of i inside and outside s , and   0 , i > 0 , and o > 0 are fixed parameters .
while the first term of the right side of ( 1 ) controls its smoothness , the surface s separates the image into distinct regions depending on their mean values .
the parameter  plays an important role in the segmentation process , as it acts like a scale parameter .
let    , s is represented as a zero level set of an implicit function  :    such that 
								 ( 2)s=={(x , y , z):(x , y , z)=0},inside(s)=={(x , y , z):(x , y , z)>0},outside(s)=={(x , y , z):(x , y , z)<0}. 
							 the level set flow associated to ( 1 ) is computed as 
								 ( 3)t=()(||i(isi)2+o(iso)2),in  ]0,[,(x , y , z,0)=0(x , y , z )  in  ,()||n=0  in   , 
							 where ( ) is the dirac function , n denotes the exterior normal to the boundary  , and 0(x , y , z ) is an initial function with the property that its zero level set corresponds to the position of the initial front . for details about using the level set method to solve ( 1 ) we refer to the original paper . 
				 
equation ( 3 ) acts only locally , on the zero level set of the function  , and the solution depends on the initial condition . in order to extend the evolution to all the level sets of  and find a global minimizer of the energy , a common approach
[ 14 , 15 ] is to set ( ) = || , leading to 
								 ( 4)t=||(ki(isi)2+o(iso)2 ) , 
							 where k =   ( ()/|()| ) is the mean curvature of every level sets of . the term k|()| performs then a regularization through a mean curvature motion .
another interesting solution has been proposed by gibou and fedkiw   who take  = 1 ( see also [ 1719 ] ) .
observing that only the zero level set of  is important and the object boundaries are represented by discontinuities in  , they ignore the regularization term and consider directly the ordinary differential equation t =   i(i  si ) + o(i  so ) .
considering the features of dataset , we suggest instead , following , to avoid any denoising and take  = 1 keeping the regularization term in ( 3 ) : 
								 ( 5)t=ki(isi)2+o(iso)2 . 
							 the first term in ( 5 )
we tested the behavior of ( 5 ) and ( 4 ) on noisy phantoms and our real data . to solve the presented pdes models we used explicit finite difference schemes as introduced in . in the next section
as previously introduced in section 2.2 , our dataset is composed by 49 volumes and represents the cell nuclei during the development of a zebrafish embryo between the sphere and the shield stage . to a first approximation ( see figure 2 ) the nuclei can be comparable to spheres of different diameter , between 6 - 7 ( internal ) and 1214 voxels ( epithelial ) . nevertheless , sometimes they are elongated and very jagged , especially before a mitosis . furthermore , the nuclei located in depth are usually incomplete .
as our aim is to segment all the nuclei at the same time , we have to segment simultaneously objects of different shapes and sizes . in the next section
we compare from this perspective the behavior of both versions of chan - vese model previously introduced .
we evaluate their performances according to different values of scale parameter  on a synthetic image and our real data .
therefore if we have to detect many objects with different size ,  should be small . if we have to detect only larger objects and to filter out smaller objects ( like points , due to the noise ) , then  has to be larger . by considering the nuclei dimension , we perform a first test on phantoms by segmenting spheres of different radii , then on real data by segmenting a region of interest selected in a volume of time lapse series .
we test both ( 4 ) and ( 5 ) using different values of . in all our computations we use i = o = 1 .
the time step is set to 0.1 if   0.5 , otherwise to 0.01 .
the segmentation error is evaluated by calculating the mean hausdorff distance ( mhd ) and the hausdorff distance ( hd )   between the segmented surfaces and a suitable gold standard . given two finite point sets , a = { a1 ,  , ap } and b = { b1 ,  , bq } , the hausdorff distance is defined as 
								 ( 6)hd(a , b)=max(h(a , b),h(b , a ) ) , 
					 where h(a , b ) = maxaa  minbb  a  b. the mhd is defined similarly , using h(a , b ) = 1/pi = 1minbb ai  b. the two formulations can be used to quantify different kinds of error : hd is related to the maximum segmentation error , whereas mhd to the total error . in figure 3
are depicted our results on noisy phantoms , using , respectively , ( figures 3(c ) and 3(d ) ) ( 4 ) and ( 5 ) .
the original data ( figures 3(a ) , 3(b ) ) are composed of spheres with radius 1 , 3 , 5 , and 7 voxels .
figures 3(c ) and 3(d )   show the original data in a cutting plane and the segmentation ( blue ) obtained using  = 0.01 , 0.1 , 0.5,1 . considering the nuclei size in our dataset
we can first observe how increasing  only bigger objects can be detected with both methods . when  = 0.01 , ( 4 ) is able to identify all the four spheres , while ( 5 ) segments also smaller regions given by noise . for 
when the scale parameter is 0.5 , only ( 5 ) still segments three regions . for  = 1 , ( 4 )
we verified that ( 4 ) is able to detect three regions only when 0.05    0.25 , while ( 5 ) can exploit a more wide range , 0.05    0.80 .
such behavior can be observed also by examining the mean and the maximum segmentation errors .
the graphs in figure 4 depict mhd and hd calculated between the surfaces of segmented spheres and corresponding surfaces of spheres created without noise .
the results are plotted only for the range of interest of . in solid line is the use of ( 4 ) . in dashed line
the sphere of the biggest radius ( 7 voxels , depicted in red ) is segmented correctly with both methods and for the whole range of scale parameter .
while using ( 5 ) errors keep about constant for all nuclei ; errors with ( 4 ) are everywhere small only for the smallest value of . the replacement of  with 1 in ( 3 ) is less sensitive to changes in size than its replacement with || .
such behavior is confirmed on real data , also if in a more restrictive useful range of values for . we selected a region of interest in a volume of our dataset .
the roi contains both epithelial and internal nuclei , and some of them are cut .
original data , gold standard , and results attained using  = 0.05 are depicted in figure 5 . for smaller values of  both methods detect regions given by noise .
we can observe that the methods described in ( 4 ) fails , because it is not able to segment all the nuclei in the subvolume .
the smallest region ( on the left in figures 5(c )   and 5(d ) ) is lost .
another small nucleus ( bottom - right in the same figures ) is not correctly segmented .
the method described by ( 5 ) is instead able to segment all the regions in the roi .
the hausdorff and the mean hausdorff distances are in this case averaged on the number of nuclei in the subvolume .
as the spacing in our data is not uniform , distances are given in m .
for  = 0.05 , 0.1 the total error is small ( 0.3 m ) and the maximum is anyhow limited ( 1.5 m ) .
the diagonal of every voxel is in fact about 1.2 m . because of robustness , we then choose to segment our time lapse series with the method described by ( 4 ) .
considering results on both phantoms and real data , we used  = 0.1 . let us point out that this analysis represents also a validation of segmentation , with very satisfactory results . 
				 
we show in this section some results obtained by segmenting the whole time lapse with the method described by ( 5 ) .
the level - set function is initialized as the difference between the 3d image i itself and its mean value .
the mean value is calculated only where i > 0 and the function thus constructed is then smoothed with few scale steps of the heat equation .
this initialization provides an initial solution close to the final one and does not depend on a previous segmentation .
actually an alternative could be to use the result of segmentation at time t as initial solution for the segmentation at time t + 1 .
but this choice does not consider that the number of nuclei changes between subsequent volumes .
indeed some nuclei divide and others are close to the bottom of the acquired volume and can go out of the imaged part of the embryo .
moreover , our solution allows the simultaneous segmentation , for example , on different processors , of all the volumes in the time lapse series .
we keep i = o = 1 , while the scale parameters have been chosen by our previous analysis and set to 0.1 .
we use in every case a time step of 0.1 and the steady state is reached after about 1500 iterations . from the isosurface corresponding to the zero level set of the final solution are then extracted the connected components , whose number provides the counting of nuclei in the processed volume .
both for extracting the isosurface zero and its connected components we make use of the open - source vtk ( visualization toolkit ) library , which provides widely used filters for image visualization and analysis .
in the first a cutting plane of last volume is compared with a cut of segmented surfaces .
the white outline well reproduces the nuclei boundaries even if they are not clearly defined , and small regions are correctly detected . in figure 7 is visualized a segmentation ( red ) together with three orthoslices of original data . observing the nuclei cut by planes it is possible to see the white contour line of original data .
as previously introduced , this number can be quantified by counting the connected components extracted from the segmented surfaces .
we show in figure 9 an example of connected components extraction on last volume of the time lapse series .
the number of nuclei detected in our dataset during the embryo development is instead plotted as graph in figure 10 . at the beginning we detected 333 nuclei ,
nevertheless we would like to point out that the number of detected cells is not strictly increasing but shows small oscillations around an exponential trend .
this behavior occurs because the cells are moving and the imaging includes only the top part of the embryo .
cells close to the bottom of the imaged volume can then be located either on the inside or on the outside . in order to estimate the ability of our method of correctly extracting the number of cells ,
we then compared the number of nuclei identified by a manual computation on the first volume of the time lapse series ( 340 nuclei ) with the number automatically detected . by our evaluation
the correct cells number is then slightly underestimated , but our error is less than 3% .
the main source of error is given by nuclei very close , if not overlapped , in original data .
this drawback could be overcome by combining signal of both nuclei and membranes , by refining the segmentation locally , or through the multiphase level - set approach .
we applied the active contour model without edges to time lapse series of confocal images representing cells nuclei in a live zebrafish embryo . as the nuclei in our dataset
are varying in shape and size , we considered two forms of the originally proposed model equation .
the maximum and total errors have been quantified as the hausdorff distance and the mean hausdorff distance between the data and a gold standard obtained through manual segmentation .
we proved that our suggestion of using  = 1 is the formulation with the best performances .
this form of the active contour model without edges has then be applied to a dataset which represents the development of a zebrafish embryo between the shield and the sphere stage . from the results of segmentation
have been afterwards extracted the connected components , whose counting allowed to quantify the number of cells .
the corresponding error has been estimated through a comparison between the number of cells detected by the algorithm and the number of cells detected by a manual computation in a selected volume . 
			
 this paper is devoted to the segmentation of cell nuclei from time lapse confocal microscopy images , taken throughout early zebrafish embryogenesis . 
 the segmentation allows to identify and quantify the number of cells in the animal model . 
 this kind of information is relevant to estimate important biological parameters such as the cell proliferation rate in time and space . 
 our approach is based on the active contour model without edges . 
 we compare two different formulations of the model equation and evaluate their performances in segmenting nuclei of different shapes and sizes . 
 qualitative and quantitative comparisons are performed on both synthetic and real data , by means of suitable gold standard . 
 the best approach is then applied on a number of time lapses for the segmentation and counting of cells during the development of a zebrafish embryo between the sphere and the shield stage . 
1. Introduction
2. Image Acquisition
3. Image Segmentation
4. Results
5. Conclusions