A reaction–diffusion based level set method for image segmentation in three dimensions

https://doi.org/10.1016/j.engappai.2020.103998Get rights and content

Abstract

The image segmentation of computed tomography data for three-dimensional biological structures remains challenging because of the limitations of existing numerical techniques and computer resources. The work represents the structures as the zero-level contour of a level set function whose value is constrained to a narrow band ranging. A cost functional composed of fitting energy for extracting the local intensity and diffusion energy for regularization is minimized within a framework of optimization. To avoid the re-initialization procedure and accelerate the convergence when updating the level set function, a reaction–diffusion technique is developed to replace the upwind algorithm by finite element analysis. Numerical examples demonstrate elegant biological structures with clear and smooth interfaces can be generated within a few iteration steps because the time step 100-fold larger than the allowable value of Courant–Friedrichs–Lewy stability condition can be applied in the proposed method.

Introduction

An efficient and clear separation method is highly desirable for three-dimensional (3D) computed tomography (CT) data because it often provides clinicians with critical diagnosis information via processed medical images. The image segmentation techniques (Haralick and Shapiro, 1985, Jin and Weng, 2019, Mandal et al., 2014), are used to partition the raw data into non-overlapped regions in accordance with their intensity or textural information. While significant progress has been made over the years in techniques to determine the tissue profiles from complex and diverse datasets, the results remain unsatisfactory because the raw data are distorted by many factors such as the sizeable inhomogeneity of tissues, artefacts of images and noise. Moreover, most computers lack resources such as powerful processors and memory to handle an enormous amount of 3D image data.

The region growing technique is one of the widely-used approaches for image segmentation in which the interfaces gradually expand from a set of given seeding points (Hojjatoleslami and Kittler, 1998, Thiran et al., 1997). In this method, the regions with similar properties such as greyscale, texture, or colour in the vicinity of the seeding points are updated in each iteration until certain predefined criteria are met. For instance, a simple yet effective hybrid seeded region growth approach is performed by assigning a seed for each distinct region (Adams and Bischof, 1994). Another well-known segmentation technique is thresholding (Mittal and Saraswat, 2018, Senthilkumaran and Vaithegi, 2016), which categorizes the elements based on the range of intensity, where the interfacial points of adjacent ranges, namely the thresholds, can be globally or locally determined. For the global criterion, the unique threshold depends on the intensity of the pre-selected points, while the multiple local thresholds rely on the neighbourhood of these points. In addition, the edge detection algorithms, originally proposed by Marr et al. (Marr and Hildreth, 1980), are also typical methods of segmenting images. Such algorithms can detect the boundaries via the convolution of intensity with various gradient-related matrixes, such as Roberts cross operator (Al-Amri et al., 2010) and Canny operator (Canny, 1987), etc. For instance, two derivative matrices can be convolved with the intensity in the horizontal and vertical directions, respectively, in order to generate the edge of a 2D figure. Sometimes, the Gaussian filter needs to be used to smooth the images before convolution (Canny, 1987).

In contrast to the above-mentioned conventional approaches, the active contour model (ACM), proposed by Kass et al. (1987), is capable of evolving initial contours to smooth and closed interfaces with sub-pixel accuracy by minimizing predefined energy functional. Under the level set framework (Osher and Sethian, 1988), most of ACMs can be categorized into region-based models (Ding et al., 2017, Min et al., 2015, Zhang et al., 2010) and edge-based ones (Li et al., 2005, Li et al., 2010, Wang et al., 2018). The region-based models utilize the image intensity to construct constraints, whereas the edge-based models apply image gradient to guide the contours toward the boundaries of desired objects. A paradigm of ACM is the well-known Chan–Vese (CV) model which determines the energy functional as the sum of the contour perimeter and fitting energy (Chan and Vese, 2001). This model can successfully segment the images into regions with distinct pixel intensity assuming there are statistically homogeneous intensities in all regions. To improve the performance of CV model in handling inhomogeneous data, by minimizing Mumford and Shah function (Mumford and Shah, 1989), Vese and Chan proposed the piecewise smooth (PS) model (Vese and Chan, 2002), where the constants, normally given as the average intensity in separated regions, are replaced by location-dependent variables to indicate the distinction of intensity in the fitting energy. It is notable that both CV and PS models are essentially region-based models because a descriptor is used to identify the regions to guide the initial contour toward real boundaries. Later, the global fitting energy in CV model is replaced by the region-scalable fitting energy which is capable of measuring the image intensities from both sides of contour. Via assigning each pixel with a homogeneous weight derived from the local greyscales distribution entropy, He et al. (2012) introduced the region-scalable fitting energy to handle intensity inhomogeneity. Moreover, to make the optimization problem well-posed, a regularization term and a penalty factor of the perimeter length are included in the objective functional (Li et al., 2008).

Even though region-based active contour models have several appealing characteristics, such as robustness to noise, insensitivity to initialization, high accuracy, and reasonable performance in dealing with intensity inhomogeneity etc., the associated enormous computational cost remains a major challenge for 3D medical images with enormous resolution. Unlike the region-based models, the edge-based active contour methods employ image gradients to terminate the contour when it approaches real boundaries. In spite of their satisfactory efficiency in segmenting most 2D images, edge-based methods demonstrate poor performance for images with massive noisy and inhomogeneous intensity. In addition, these methods often become less efficient for 3D medical data because of complicated images characterized with blur interfaces, boundary leakage and difficulties in controlling the motion of contour.

It is noteworthy that most previous studies in the image processing adopt the upwind finite difference algorithm (Huang and Wang, 2008, Wang et al., 2000) to update the level set function when solving Hamilton–Jacobi equation. However, the major shortcomings of such algorithms such as the small time step owing to Courant–Friedrichs–Lewy (CFL) stability condition (De Moura and Kubrusly, 2013) and frequent re-initialization of level set function for the sake of regularity (Osher and Fedkiw, 2006) have not been solved well. Thus, the development of an efficient algorithm to process an enormous amount of voxel-based data obtained from high-resolution CT images of sophisticated biological tissues is highly desirable.

Fortunately, some methods have been proposed to increase the time step size and eliminate the re-initialization procedure. For instance, based on the operator splitting techniques (Weickert et al., 1998), the computation is only performed on a narrow band (Ye et al., 2012) around the given surface dataset so that the limitation of the CFL condition is removed (Li and Kim, 2015). Additionally, by applying the linearly stabilized splitting scheme (Eyre, 1998) with the Fourier-spectral method (Sheng et al., 2010), a larger time step is allowed in an image reconstruction technique (Li et al., 2019). Nevertheless, expensive re-initialization cannot be completely abandoned for the sake of regularity in these approaches. In conventional level set methods, re-initialization has been extensively employed as a remedy for avoiding numerical errors and maintaining the stability of contour evolution, but this procedure requires heavy computational time. For overcoming this difficulty, Li et al. (2010) introduced a distance regularizing term in the energy functional to intrinsically retain the regularity of level set function during the level set evolution. Then, by combining the reaction–diffusion technique with the level set method for image processing, a new level set evolution approach becomes completely free of the costly re-initialization procedure (Zhang et al., 2013). Such an approach adopts a two-step splitting method to solve the energy functional iteratively, which first iterates the level set evolution equation and then solves the diffusion equation. Instead of re-initialization, the level set function is regularized by a diffusion equation from the reaction–diffusion model. But the CFL condition associated with the upwind algorithm is still a difficulty in these methods. Thereafter, by introducing an additive operator splitting algorithm and incorporating a signed distance preserving term for regularizing the evolution, a distance preserving level set method is proposed to avoid the re-initialization and the upwind algorithm (Estellers et al., 2012).

This study develops a reaction–diffusion based level set method to expedite the update of the level set function, which involves a reaction process to describe local interaction between intensities and diffusion to allow the intensities to be homogenized in a spatial context. The idea of the reaction–diffusion has been successfully introduced to structural topology optimization (Choi et al., 2011, Yamada et al., 2010) to enable the convergence of objective function within a few iteration steps. A more attractive feature of using the reaction–diffusion method resides in the indispensability of the re-initialization procedure by intercepting the level set function at a range of −1 and 1 and regularizing the level set function with a diffusion term. In this method, the regularity of the contour can be fully maintained by the diffusion term, thus it is unnecessary to include the term of perimeter length to further improve the regularity. Moreover, the paper uses a finite element method to update the level set function, which allows the time step size many times greater than the criterion of the CFL condition.

The rest of the paper is organized as follows. In Section 2, the optimization methodology and the proposed updating technique are reviewed. Several numerical examples with and without normalization and the corresponding discussion are presented in Section 3. In Section 4, the conclusion is drawn.

Section snippets

Level set method

In this paper, the interfaces Γ of images are represented by the zero-level contour of a high-dimensional level set functional in Lipschitz space. Mathematically, it is represented as: φ(x)>0xΩ1φ(x)=0xΓφ(x)<0xΩ2where x denotes the spatial variable, and the level set function φ separates the domain Ω=Ω1 Ω2 Γ into void region Ω1 and solid region Ω2 as per its value. Note that φ(x) is replaced by φ in the following context for simplicity. The evolution of the boundaries with respect

Numerical results and discussion

In this section, the reaction–diffusion based level set method is employed to segment a woodpecker’s CT images with a voxel as high as 120 × 350 × 153 for assessing the validity and applicability of the proposed method. The raw data is obtained from a positron emission tomography-computed tomography (PET-CT) scanner with a spatial resolution lower than 1.3 mm and provided by a high-tech molecular imaging device enterprise (Pingseng-Healthcare, 2012). The image slice of tomography data on the

Conclusion

This paper proposed a reaction–diffusion based level set method for segmenting 3D biological structures from high-resolution raw CT data in this paper. It provides a novel solution to the difficulties in massive existing separation methods, such as frequent re-initialization and small time step. As the diffusion energy is included in the cost function, the regularity of the level set function and the smoothness of the contour can be intrinsically maintained and therefore the re-initialization

CRediT authorship contribution statement

Zhe Zhang: Conceptualization, Methodology, Software, Writing - review & editing. Yi Min Xie: Formal analysis, Writing - review & editing, Supervision. Qing Li: Writing - review & editing. Shiwei Zhou: Conceptualization, Methodology, Software, Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

The work was supported by the Natural Science Foundation of Shanghai, China (No. 18ZR1440900).

References (43)

  • MinH. et al.

    An intensity-texture model based level set method for image segmentation

    Pattern Recognit.

    (2015)
  • MittalH. et al.

    An optimum multi-level image thresholding segmentation using non-local means 2D histogram and exponential Kbest gravitational search algorithm

    Eng. Appl. Artif. Intell.

    (2018)
  • OsherS. et al.

    Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations

    J. Comput. Phys.

    (1988)
  • SethianJ.

    Evolution, implementation, and application of level set and fast marching methods for advancing fronts

    J. Comput. Phys.

    (2001)
  • ThiranJ. et al.

    A queue-based region growing algorithm for accurate segmentation of multi-dimensional digital images

    Signal Process.

    (1997)
  • WangL. et al.

    Active contours driven by edge entropy fitting energy for image segmentation

    Signal Process.

    (2018)
  • YamadaT. et al.

    A topology optimization method based on the level set method incorporating a fictitious interface energy

    Comput. Methods Appl. Mech. Eng.

    (2010)
  • ZhangK. et al.

    Active contours driven by local image fitting energy

    Pattern Recognit.

    (2010)
  • AdamsR. et al.

    Seeded region growing

    IEEE Trans. Pattern Anal. Mach. Intell.

    (1994)
  • Al-AmriS. et al.

    Image segmentation by using edge detection

    Int. J. Comput. Sci. Eng.

    (2010)
  • ChanT. et al.

    Active contours without edges

    IEEE Trans. Image Process.

    (2001)
  • Cited by (9)

    • Study on the effect of extreme learning machine and its variants in differentiating Alzheimer conditions from selective regions of brain MR images

      2022, Expert Systems with Applications
      Citation Excerpt :

      Active contour models have been employed in the segmentation of these brain regions (Anandh, Sujatha, & Ramakrishnan, 2014; Ramaniharan et al., 2016). Level set methods that utilize curves or surfaces and partial differential equations are known to be effective for MR image segmentation (Zhang, Xie, Li, & Zhou, 2020). In order to characterize the segmented brain regions, different feature extraction techniques are employed.

    • Learning finite difference methods for reaction-diffusion type equations with FCNN[Formula presented]

      2022, Computers and Mathematics with Applications
      Citation Excerpt :

      Among various PDEs representing natural phenomena, we consider reaction-diffusion type equations. The reaction-diffusion model has been applied and used in various fields such as biology [3–6], chemistry [7–9], image segmentation [10–12], image inpainting [13–15], medicine [16–18]. In this paper, we focus on second order reaction-diffusion type equations, including the heat, Fisher's, Allen–Cahn (AC) equation [19], and reaction-diffusion equations with trigonometric function terms.

    • A reaction diffusion-based B-spline level set (RDBLS) method for structural topology optimization

      2022, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      In addition, the geometrical complexity of the optimized structure could be flexibly controlled by incorporating fictitious interface energy. Based on these advantages, the reaction diffusion-based level set method (RDLS) has been applied in multi-material structures [49], image processing [50], and further developed within the body-fitted mesh. Specifically, the body-fitted mesh was utilized in RDLS for structural optimization to obtain ultra-smooth boundaries [24].

    • Body-fitted bi-directional evolutionary structural optimization using nonlinear diffusion regularization

      2022, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      Besides, the elegant geometric representation can be obtained directly without post-processing or material interpolation. Compared to our previous work, we found the solid–void boundaries based on the relative density distribution instead of the level set function [45]. Let us assume that Fig. 3(a) represents the density field over a rectangular design domain.

    • A reaction diffusion-based level set method using body-fitted mesh for structural topology optimization

      2021, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      Unlike the traditional upwind scheme [57–59], the restraint to the time step due to the Courant–Friedrichs–Lewy (CFL) stability condition is invalid in the RDLS method [60,61]. Thus, the large time step can be used to reduce the number of iterations [62]. The fictitious interface energy term, together with the interception of the level set function in {-1,1} regularizes the level set function, making the level set function re-initialization unnecessary [53].

    View all citing articles on Scopus
    View full text