4D Segmentation Algorithm with application to 3D+time Image Segmentation

11/21/2021
by   Markjoe Olunna Uba, et al.
0

In this paper, we introduce and study a novel segmentation method for 4D images based on surface evolution governed by a nonlinear partial differential equation, the generalized subjective surface equation. The new method uses 4D digital image information and information from a thresholded 4D image in a local neighborhood. Thus, the 4D image segmentation is accomplished by defining the edge detector function's input as the weighted sum of the norm of gradients of presmoothed 4D image and norm of presmoothed thresholded 4D image in a local neighborhood. Additionally, we design and study a numerical method based on the finite volume approach for solving the new model. The reduced diamond cell approach is used for approximating the gradient of the solution. We use a semi-implicit finite volume scheme for the numerical discretization and show that our numerical scheme is unconditionally stable. The new 4D method was tested on artificial data and applied to real data representing 3D+time microscopy images of cell nuclei within the zebrafish pectoral fin and hind-brain. In a real application, processing 3D+time microscopy images amounts to solving a linear system with several billion unknowns and requires over 1000 GB of memory; thus, it may not be possible to process these images on a serial machine without parallel implementation utilizing the MPI. Consequently, we develop and present in the paper OpenMP and MPI parallel implementation of designed algorithms. Finally, we include cell tracking results to show how our new method serves as a basis for finding trajectories of cells during embryogenesis.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 7

page 10

page 11

page 12

page 32

page 33

page 35

page 36

06/22/2020

A Proximal-Gradient Algorithm for Crystal Surface Evolution

As a counterpoint to recent numerical methods for crystal surface evolut...
10/03/2007

Colour image segmentation by the vector-valued Allen-Cahn phase-field model: a multigrid solution

We propose a new method for the numerical solution of a PDE-driven model...
04/19/2021

Two-phase image segmentation by the Allen-Cahn equation and a nonlocal edge detection operator

Based on a nonlocal Laplacian operator, a novel edge detection method of...
07/17/2020

Anisotropic Mesh Adaptation for Image Segmentation Based on Mumford-Shah Functional

As the resolution of digital images increase significantly, the processi...
06/17/2020

A parallel hybrid implementation of the 2D acoustic wave equation

In this paper, we propose a hybrid parallel programming approach for a n...
08/25/2017

Linear Differential Constraints for Photo-polarimetric Height Estimation

In this paper we present a differential approach to photo-polarimetric s...
11/01/2010

Volume-Enclosing Surface Extraction

In this paper we present a new method, which allows for the construction...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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 literature

segm5 ; 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

(1)

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:

(2)

Equation (1) is accompanied by Dirichlet boundary conditions

(3)

or Neumann boundary conditions

(4)

where is unit normal to , and with the initial condition

(5)

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

(6)

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.

Figure 1: First column shows a solid sphere and sphere with six symmetric holes, whereas the second column shows results of the segmentation after the application of the model (1).
Figure 2: First column of this figure shows the 3D microscopy image of cell membrane to be segmented, whereas the second column shows the result after application of the model (1); the cell membrane of interest and its corresponding result after segmentation (which are depicted in black) is located approximately at the center of each picture.
Figure 3: First column of this figure shows the 3D microscopy image of cell nuclei together with approximate cell centers (these are points marked in red color), whereas the second column shows the 3D image of microscopy image of cell nuclei after application of (1) to the cells with centers marked in red color.
Figure 4: First column of this figure shows the 3D microscopy images of cell nuclei to be segmented, whereas the second column shows the result after application of the model (1); the cell nuclei of interest and their corresponding results after segmentation (which are depicted in black) are located approximately at the center of each picture.

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.

Figure 5: This figure shows the D microscopy image of mouse embryo cell nuclei segmentation result using classical SUBSURF approach.
Figure 6: First column of this figure shows the 3D microscopy images of cell nuclei to be segmented, whereas the second column shows the result of segmentation using the classical SUBSURF approach.
Figure 7: First row of this figure shows the 3D microscopy images of cell membranes to be segmented, while the second row shows the result of segmentation using the classical SUBSURF approach.

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 57 using this approach, we obtained the results presented in Figures 810, which is clearly much more accurate compared to the results, shown in Figures 57, obtained using the classical method.

Figure 8: This figure shows the D microscopy image of mouse embryo cell nuclei.
Figure 9: In this figure, first row shows the D microscopy images of four different cell nuclei; second row shows segmentation result using thresholded image intensity information and original image intensity information in 3D case of (7). That is, the result after application of the 3D case of (7) with and .
Figure 10: First column of this figure shows the D microscopy images of cell membrane to be segmented. The second, third, and fourth columns show the result after the application of (7) in 3D with and , and , and and , respectively.

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

where .

Our new method is based on solution of the following generalized SUBSURF equation

(7)

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

(8)

and with the initial condition

(9)

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

(10)

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

(11)

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

(12)

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

(13)

Equation (13) can be rewritten as

(14)

which can be written in the form of system of equations

(15)

, , , and . The 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

(16)

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)

(17)

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

(18)

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:

(19)

, , , and . The coefficients of (19) are computed using instead of .

Finally, we note that if , then , , and , and the equation (19) simplifies to

(20)

Since in our implementation we used common , in the sequel, we deal with the properties of the system (20). All derived properties are simply adopted also to the system (19).

The solution of the linear system given by equations (20) is obtained by the successive overrelaxation (SOR) method as follows:

(21)

where