In image processing and computer vision, it is well known that image segmentation is one of the fundamental and most studied problems. Consequently, there are several approaches to image segmentation in literaturesegm5 ; segm2 ; segm3 ; segm4 ; segm1 ; segm7 ; segm8 ; sethianb ; segm6 . However, in this paper, we will study and use the level set method for image segmentation. The idea of following the interfaces in immiscible fluid flow by the Eulerian approach was proposed by Dervieux and Thomasset levelset1 ; levelset2 , and Osher and Sethian oshersethian then introduced the general level set methodology. The level set method is the evolution of curves and surfaces by means of a dynamical embedding function (often referred to as the level set function). In other words, the level set method is the implicit representation of boundaries, curves (or surfaces). Furthermore, in the level set formulation, the subjective surface (SUBSURF) method for image segmentation will be the focus of this paper. This method, in the context of image processing, was introduced in sartietal ; subsurf , studied and applied in several biomedical research cmssga ; faure ; app1 ; balasz ; kmnapmr ; kmikula ; mprs ; app4 ; app2 ; sartietal ; subsurf ; app3 . SUBSURF segmentation method is based on the idea of evolution of segmentation function governed by a nonlinear diffusion equation cmssga ; sartietal ; subsurf ; zcrmbmpa which can be understood as an advection-diffusion model mprs . Hence, a segmentation seed (the starting point that determines the approximate position of an object in the image) is usually needed to segment an object. Then an initial segmentation function is constructed with reference to the segmentation seed. Finally, this segmentation function is allowed to evolve to the final state following the SUBSURF model. Ideally, the evolution process ends up with a function whose isosurfaces all have the object’s shape that is intended to be segmented. At each time step in the evolution process, the values of the segmentation function are locally rescaled to interval . After the last time step in the evolution, the contour is chosen as the approximate boundary of the segmented object.
In this paper, we introduce a generalization of the classical SUBSURF model. This generalization is due to the fact that in real applications where the object intended to be segmented possesses internal structures or edges, it is usually challenging to obtain optimal results using the classical SUBSURF segmentation approach. The reason is that this approach works with edge information throughout the segmentation process. Hence, edges within the internal structures in an object of interest are also respected during segmentation. To overcome the effect of the internal structures or edges, we introduced thresholding of image intensity values within a ball of appropriate radius around the object center. This local thresholding serves to eliminate the internal structures or edges. Additionally, we combined the information obtained from thresholding and original image intensities to get an accurate final segmentation result. Unlike the works in kmikula ; mprs ; app4 ; msarti ; app2 ; app22 , at each segmentation step, the values of the segmentation function are rescaled to interval . Finally, considering that the computation task involving 4D microscopy images requires solving a linear system with several billion unknowns, we developed an efficient parallel implementation of the method for the generalized model in the case of 4D images.
The remaining part of this work is organized as follows: Section provides an overview of the SUBSURF method and some numerical examples of the application of the SUBSURF segmentation method to artificial and real data. Section contains the details of the proposed generalization of the classical method. It contains the studied model combining thresholded image information and original image data, the numerical scheme for solving the model, the stability of the numerical scheme, experimental order of convergence, an overview of the computer implementation, and parallel implementation using OpenMP and MPI. Furthermore, the application of the new segmentation algorithm to D image segmentation was presented in Section . Finally, cell tracking based on the result of the new segmentation algorithm was presented in Section .
2 The SUBSURF segmentation method
The SUBSURF segmentation method is based on the idea of evolution of segmentation function, which is governed by an advection-diffusion model. The SUBSURF model was proposed by Sarti, Malladi, and Sethian subsurf and is given by
where is the unknown segmentation function, , is image to be segmented, and is a smoothing kernel. is a smooth nonincreasing function such that and as . The function was introduced by Perona and Malik in pm and has the following typical example:
Equation (1) is accompanied by Dirichlet boundary conditions
or Neumann boundary conditions
where is unit normal to , and with the initial condition
Dirichlet boundary condition is used in the image segmentation, and without loss of generality may be assumed. The zero Neumann boundary
conditions are used, e.g., in morphological image smoothing (see,
e.g., cmssga and references therein).
In computations involving model (1), one possible choice of initial segmentation function, , is the following function
where corresponds to the focus point, and gives a maximum of whose peak is centered in a “focus point” inside the segmented object app2 . This function can also be considered only at a circle with center and radius . Outside this circle, the value of is taken to be . If (1) is accompanied by zero Dirichlet boundary condition, then finally is subtracted from such peak-like profile. In the case of small objects, a smaller choice of value for can be used to speed up computations. Also, is another possible choice of initial segmentation function.
We now present some numerical and illustrative examples to demonstrate the application of the SUBSURF model to image segmentation. The model (1) for surface reconstruction was tested on artificial examples and applied to real data representing 3D microscopy images of cell nuclei. For the first numerical experiment, a sphere and a sphere with holes were generated, as can be seen in the first column of Figure 1. Afterward, model (1) was used to reconstruct the shapes of these spheres, and the results obtained after reconstruction are shown in the second column of Figure 1. In the second experiment, the model (1) was applied to real data representing 3D microscopy images of cell nuclei and membrane. In the first column of Figures 2, 3, and 4, 3D image of these microscopy images are shown, whereas, in the second column of Figures 2, 3, and 4 the results obtained after application of model (1) are shown and colored in black.
3 4D Image Segmentation Algorithm
In the previous section, we presented the SUBSURF segmentation method and some of its successful applications. However, in many real applications where the object to be segmented has internal structures or edges, it is usually challenging to obtain optimal results using the SUBSURF segmentation approach (see, e.g., Figures 5, 6, and 7). The reason is that this approach works with edge information throughout the segmentation process. Hence, edges within the internal structures in an object of interest are also respected during segmentation. For instance, Figures 5, 6, and 7 show the results of segmentation using the classical SUBSURF approach for D images of zebrafish cell nuclei, cell membrane, and mouse embryo cell nuclei, respectively. The results of segmentation in these figures are colored in black, and the shapes of interest are located approximately at the center of each image.
In this section, to overcome the effect of the internal structures or edges, we introduce the thresholding of values within a ball of appropriate radius around the object center. This local thresholding serves to eliminate the internal structures or edges. Finally, we combine the information given by thresholding and original image intensities to get a segmentation result. For example, in D segmentation (see also ubaetal ; ubaetalm ) of nuclei and membranes images shown in Figures 5–7 using this approach, we obtained the results presented in Figures 8–10, which is clearly much more accurate compared to the results, shown in Figures 5–7, obtained using the classical method.
In 4dpaper , the first numerical scheme for D image segmentation based on a generalized SUBSURF model mprs was introduced. In this paper, we introduce and study a D numerical scheme for the solution of the new D segmentation model (7) presented in the next section. This numerical method is based on the finite volume approach. The stability of the numerical scheme is also presented. Additionally, we introduced the local rescaling of the values of the segmentation function at each segmentation step to the interval . Furthermore, the application of the studied model to D and D image segmentation is studied. OpenMP and MPI parallelizations are employed for efficient computer implementation. For example, in the real application presented in the fifth and sixth experiments of Section 4, to process the 3D+time microscopy images, with dimensions , 1012 GB of memory was needed. This clearly shows that it may not be possible to process these images on a serial machine without parallel implementation utilizing the MPI.
3.1 Mathematical model
Let , be the intensity function of a 4D image where , , and represents the final real video time. For each real video time , let denote the set of “centers,” where and represent the center and the total number of centers, respectively, at time . Furthermore, for each and , let , , where is the ball at time with radius centered at , a given point (“center”) inside the object to be segmented. Then the threshold value (which is used for local thresholding) may be chosen as , and the ball radius may be chosen with respect to the approximate size of the object to be segmented. So, the idea (of local thresholding) is to set all intensity values in the local neighborhood of center to if they are above and otherwise. Hence, the D image intensity of the thresholded image within a ball of radius is defined by
Our new method is based on solution of the following generalized SUBSURF equation
where is the unknown segmentation function, , is a “segmentation time,” , is the function introduced by Perona and Malik in pm (see, e.g., equation (2)), determine the influence of information obtained from thresholded and original image intensities, is the smoothing kernel, e.g., the Gaussian function, and are presmoothing of and by convolution with , respectively. It is well known, see e.g., witkin ; koend , that the convolution with Gaussian function is equivalent to solving the linear heat equation (LHE) obtained by setting in equation (7), with . Thus, the presmoothing of and by convolution with Gaussian function is realized by solving the LHE. Equation (7) is accompanied by Dirichlet boundary conditions
and with the initial condition
Without loss of generality, is assumed. The initial condition in D is constructed using equation (6), where .
3.2 Numerical discretization
In this section, the details of time and space discretization is presented.
3.2.1 Time discretization
For time discretization of (7), semi-implicit approach which guarantees unconditional stability is used. Suppose that the (7)–(5) is solved in time interval and equal number of time steps. If denotes the time step, then the time discretization of (7) is given by
where is the regularization factor (Evans - Spruck evans ), is given initial segmentation function, and , is the solution of the model in segmentation time step .
3.2.2 Space discretization
For space discretization, we start with introduction of some notations which will be used subsequently. We have adopted similar notations as those used in mprs and kmikula . Let denote finite volume mesh containing the doxels of D image, while , , , , denote each finite volume. For each , let , , , be the size of the volumes in , , , direction. Let the volume of and its barycenter be denoted by and , respectively. Let the approximate value of in be denoted by . For every , we have the following definitions:
For each , denote the line connecting the center of and the center of its neighbor by and its length . We denote the sides, its measure, and normal in coordinate directions of the finite volume by , , and , respectively. Furthermore, for each , is defined by
The approximate value of in ,
with , is denoted by ; the time index is omitted because only the values from the time level will be needed at these points.
With these notations, integration of (10) over finite volume yields
Let the average value of in finite volume be denoted by . If we consider the fact that and are asummed to be piecewise constant over the finite volume mesh on the left hand side of (11) and using the divergence theorem on the right hand side of (11), we obtain
where beneath the summation sign, we used just instead of to simplify notation. If we approximate normal derivative by and define and to be the average of and on then (12) reduces to
Equation (13) can be rewritten as
which can be written in the form of system of equations
, , , and .
average values , , and either in doxels or on doxel sides
are determined using the reduced diamond cell strategy (see kmikula ) adapted to D.
In the sense of this reduced diamond cell approach, the approximate values of are obtained in the points . These values are defined for each by
The components of the averaged gradient on , , are approximated by D diamond cell approach which use the values given above (see also kmikula ). Additionally, approximation of the gradient on the face is denoted by . This implies that
If the same approach for computation of gradients of image intensities is used, then the following approximation of in (15) is obtained
where , is the value of in doxel , and , is the value of in doxel . Finally, incorporating regularization, we obtain the following remaining terms in (15)
Equations (15) accompanied by the zero Dirichlet boundary condition, represent a linear system of equations which can be solved efficiently, e.g., using the Successive Overrelaxation (SOR) method.
Remark 1 (Local rescaling).
For each segmentation time step , and , let denote the set of centers, where and represent the center and the total number of centers, respectively for each . Let be a ball with radius and center ( is a given point inside the object to be segmented), and . Then the locally rescaled version of given by (15) within the ball is obtained by the following relation
Consequently, we have that for each time step , rescaled version and it is used instead of in (15) in the next segmentation time step. So in each segmentation step we solve the following system of equations representing the final formulation of our method:
, , , and . The coefficients of (19) are computed using instead of .
Finally, we note that if , then , , and , and the equation (19) simplifies to
The solution of the linear system given by equations (20) is obtained by the successive overrelaxation (SOR) method as follows: