# A Markov Random Field and Adaptive Regularization Embedded Level Set Segmentation Model Solving by Graph Cuts

Authors: Dengwei Wang
Journal: Journal of Electrical and Computer Engineering
 (2019)

---

## Abstract

This paper presents a novel Markov random field (MRF) and adaptive regularization embedded level set model for robust image segmentation and uses graph cuts optimization to numerically solve it. Firstly, a special MRF-based energy term in the form of level set formulation is constructed for strong local neighborhood modeling. Secondly, a regularization constraint with adaptive properties is imposed onto the proposed model with the following purposes: reduce the influence of noise, force the power exponent of the regularization process to change adaptively with image coordinates, and ensure the active contour does not pass through the weak object boundaries. Thirdly, graph cuts optimization is used to implement the numerical solution of the proposed model to obtain extremely fast convergence performance. The extensive and promising experimental results on wide variety of images demonstrate the excellent performance of the proposed method in both segmentation accuracy and convergence rate.


---

## Body

## 1. Introduction

Image segmentation is the technology and process of dividing an image into several specific regions with unique properties. In the field of computer vision, it is a key step from image processing to image analysis. Researchers in image processing direction have proposed a large number of algorithms for image segmentation; however, none of them can solve the image segmentation problems of all modalities; therefore, the image segmentation direction has been developing. Among the family of image segmentation algorithms, there is a kind of method that occupies an absolute advantage, that is the geometric active contour models based on level set representation; such methods have a very rigorous mathematical foundation and are capable of handling the topological changes of the contour freely, which are difficult to be solved by the parametric active contour models, and thus it can simultaneously segment multiple targets and get a smooth and closed target contour. Without loss of generality, we can roughly classify the existing LSMs into the following two categories: the edge-based models [1, 2] and the region-based models [3–8].The edge-based models use the gradient information of the image to construct the driving force required for the evolution process. The geodesic active contour (GAC) model [1] proposed by Caselles et al. is a typical representative, and it is also one of the most successful and popular edge-based models; the driving force of its segmentation is derived from the intrinsic geometric measure of the input image. In order to completely eliminate the time-consuming reinitialization operation, Li et al. [2] presented a new level set regularization (LSR) energy term, which is used to penalize the deviation of the current level set function from the signed distance function, and embedded it into a variational level set model with the edge-based stopping energy term.The region-based modes construct the driving force needed for the evolution process based on the regional statistical information of the image. The popular Chan–Vese (CV) model [3] which was proposed on the assumption that image intensities is statistically homogeneous (roughly a constant) in each region. By transforming the problem of curve evolution into a simple integer operation between two linked lists, Shi and Karl [4] constructed a fast two-cycle (FTC) algorithm for the approximation of level-set-based curve evolution. In order to solve the problem of inhomogeneous image segmentation, Li et al. [5] proposed a variational level set model based on local binary fitting (LBF) energy. Lankton and Tannenbaum [6] proposed a natural framework called localized active contours (LAC), which allows any region-based level set segmentation models to be reformulated in a local way. Ding et al. [7] designed a local prefitting energy- (LPFE-) based active contours model for fast image segmentation, which achieves relatively fast segmentation performance and relatively robust initialization performance. In addition, relying on reciprocal cross-entropy theory, Ni and Wu [8] introduced an active contours model for image segmentation applications based on the novel fitting term (NFT).The aforementioned two broad categories of methods can obtain ideal segmentation results in the high-quality image segmentation applications; however, the segmentation performance of these methods drops sharply when the noise interference components in the test image are obvious. The reason for this is that the existing methods usually assume that the pixels in each region are independent of each other when constructing the energy functionals; such a hypothetical logic makes the curve evolution process very sensitive to noise information. In addition, most of them use traditional numerical schemes of finite difference type to discretize the segmentation models. However, the internal logic and execution mechanism of the traditional numerical schemes of finite difference type determine its time consumption is very serious, thus this type of numerical strategy is powerless for image segmentation applications with high real-time requirements.To solve the above problems, this paper presents a novel MRF [9] and adaptive regularization embedded level set model for robust image segmentation and uses the graph cuts to numerically solve it. Firstly, a special MRF-based energy term in the form of level set formulation is constructed for strong local neighborhood modeling. With the help of MRF modeling, adjacent pixel groups are preferentially assigned to the same region, which can be guaranteed even in the presence of noise interference, thus it has natural robustness to noise interference. Secondly, a regularization constraint with adaptive properties is imposed onto the proposed model with the purpose of reducing the influence of noise, ensuring that the active contour does not pass through the weak object boundaries. The adaptive regularization here is the extension and optimization of the method proposed by Zhou and Mu [10]. After extension, the proposed model overcomes the problem that Zhou and Mu’s method can only take a fixed constant exponent and cannot change adaptively with the pixel position. Thirdly, the graph cuts optimization [11] is used to implement the numerical solution of the proposed model. Under its help, our model obtains very fast convergence performance, and this can be confirmed by the quantitative data in the experimental section.The remainder of this paper is organized as follows: Section2 is a brief description of the related background knowledge. Section 3 presents the proposed model and its corresponding implementation strategies. Section 4 validates the proposed model by extensive experiments and discussions on a wide variety of images. Last, conclusions are drawn in Section 5.


## 2. Background

The LSM was originally proposed by Osher and Sethian [12] in 1988 to solve the change process of the flame shape that follows the thermodynamic equation. This is due to the high dynamism of the flames and the variability of the topological structure; if parametric curves or surfaces are used to describe this change of flame, it will inevitably encounter great difficulties. However, this problem can be solved well by introducing the level set function framework.The core idea of level set image processing is to represent an active contour as the zero level set of an implicit function (called level set function) which is one dimension higher than it. Figure1 shows an example of a level set function and its corresponding zero level set.Figure 1
Level set function and its corresponding zero level set.When dealing with the evolution of the planar curves, the LSM is not trying to track the location of the evolving curve, but follows a certain law (usually a gradient descent flow equation); in the two-dimensional Cartesian coordinate system, the LSF is updated continuously so as to achieve the goal of evolving the closed curve which is hidden in the LSF. The biggest advantage of this curve evolution style is that the LSF remains a valid function even if the topological change (splitting or merging) occurs in the closed curve which is hidden in the LSF. Figure2 shows an example of the LSE; the first row represents the three states in the LSE process, and the second is the zero level set curve corresponding to the first row. From this diagram, we can find that the topological structure changes of the evolving curve can be well handled (the topological change here is splitting type) by using the level set expression.Figure 2
An example of LSE. (a) Level set function and zero horizontal plane. (b) The zero level set corresponding to the level set function.Below we give a mathematical description of the level set function: the LSM represents the moving curve (or active contour)Ct by the zero level set of a Lipschitz function ϕt,x,y, i.e., Ct=x,yϕt,x,y=0. The evolution equation of the level set function ϕ can be written as the following general form:(1)∂ϕ∂t+F∇ϕ=0,where “∇” is the gradient operator and ϕ represents the iteration velocity of the level set evolution process. For image segmentation applications, the speed function F is usually determined jointly by the image data and the level set function ϕ.


## 3. The Proposed Model

In this section, we will present the segmentation model proposed in this paper, which uses MRF and adaptive regularization as its kernel of thought. Similar to most popular LSMs, here we represent the energy functional framework of the proposed segmentation model as follows:(2)E=EREG+EMRF,where EREG is the energy term used to implement the regularization operation and EMRF is the energy term that models the relationship between the pixel and its neighborhood under the Markov random field framework. At the macro role level, EREG is equivalent to the internal energy and its function is to modify the properties of level set function itself. EMRF is equivalent to the external energy, and its role is to promote the level set function toward the direction of the foreground.Below, we will first construct the energy subitems that the proposed model relies on. Finally, give a detailed solution process based on the idea of graph cuts optimization.


### 3.1. Energy Terms

#### 3.1.1. Regularization Term

In order to enhance the smoothness of the level set evolution process, we introduce an adaptive regularization constraint that is directly applied to the zero level set, the effect of which is to reduce the influence of noise and ensure that the active contour does not pass through the weak object boundaries.In order to solve the problem of regularization of zero level curve, Zhou and Mu [10] proposed an efficient regularization scheme based on weighted Dirichlet integral. Although this method has achieved good results, its constant exponent gives it the following problems: cannot effectively reflect the local characteristics of the image; therefore, the regularization scheme cannot let the exponent automatically adapt to the image data. In order to enhance the adaptability of the constraint to the local image information, we propose the following weighted Dirichlet integral regularization with variable exponent:(3)EREGϕ=∬Ωδϕ∇ϕf∇Gσ∗If∇Gσ∗Idxdy,where “Gσ∗I” is the convolution operation between the image I and the Gaussian filter kernel Gσ (whose standard deviation is σ), “∇” is a gradient operator, and the exponent fn is a strictly monotone function and has the following limiting properties: (i) limn⟶0fn=1; (ii) limn⟶+∞fn=1.5. By matching these two features, we can construct the following simple function:(4)fn=1.5−1n+2.Below, we briefly analyze the macroscopic properties of the functionEREG−ZEROϕ:(i)
Whenf∇Gσ∗I equals to 1, RLSCϕ will degenerate to the following simple form: RREGϕ=∬Ωδϕ∇ϕdxdy. By comparing it with the C-V model, we find that this degenerate form is the curve length constraint used in the C-V model, which calculates the length of the zero curve of the level set function ϕ.(ii)
Whenf∇Gσ∗I is a constant F greater than 1, EREGϕ will be simplified as a geometric constraint for zero level curve used in the literature [10].For most images, the intensity value of the same foreground area is not strictly uniform, that is to say, different local areas within the foreground range have different local characteristics. Under this generality, if the above exponentf⋅ is set as a constant, the local characteristics of the image cannot be correctly reflected; in other words, the exponent f⋅ cannot automatically match the image data. Therefore, it is very unreasonable to set the exponent f⋅ constant. In view of this, we take a different approach. In our framework, the exponent f⋅ is directly related to the local intensity information of the image. Our approach ensures that weak regularization (f∇Gσ∗I≈1) is performed within regions with almost constant intensity values (the gradient value of the image is almost zero); as a result, it can effectively avoid the problem of the disappearance of weak boundaries. In contrast, in other areas, it will take a strong regularization operation to force the removal of false contours. Therefore, the existence of variable exponent f⋅ makes our regularization scheme adaptively choose between weak regularization and strong regularization. In the experimental results section, we will through a simple example to intuitively feel the effect of the variable exponent strategy.


#### 3.1.2. MRF-Based External Energy Term

The traditional level set segmentation models ignore the coupling relationship between neighboring pixels, so the evolution curve often converges to the wrong noise position. In order to improve the overall antinoise performance of the level set segmentation model, we introduce the idea of MRF to model the coupling relationship between the pixels and their neighborhoods.Firstly, we represent the square neighborhood required by the MRF idea asℕwm, whose size is 2w+1×2w+1 and m refers to the component of the random variable M. In order to match the requirements of MRF, we need to first construct a coupled random field whose purpose is to design our variable and symbolic system, here we represent the coupled random field as CRF=I,C, where I=Im,m∈M represents the field of pixel intensity, and C=Cm,m∈M represents the label random field (also known as characteristic random field), which is used to distinguish the foreground (F) from the background (B).In addition, the relevant theoretical knowledge of MFR tells us that the energy functional based on MRF finds the optimal category (foreground or background) label by using the neighborhood information of each pixel. Usually, we can achieve this goal by maximizing a posteriori (MAP) probability distribution function. With the help of Bayes theory, we can expressPCI as follows:(5)PCI=PIC⋅PCPI,where PCI, PC, PIC, and PI represent segmentation probability, priori segmentation probability, conditional segmentation probability, and image probability, respectively. When the image data are given, the value of PI is a constant; therefore, we can delete it from the above expression, and the deletion will not affect the proportional relationship contained in the equation. Therefore, we have the following further proportional relationship:(6)PCI∝PIC⋅PC.The existing level set segmentation models generally assume that the pixels in the foreground and background regions are independent from each other. However, such assumption is not reasonable and directly leads to their segmentation results vulnerable to noise interference. In order to overcome this shortcoming, we introduce the MRF theory to model the two components ofPIC and PC in the above expression, that is to say, when we construct the probability expression, we incorporate the neighborhood relationships into our consideration system. After the introduction of MRF, the values of PIC and PC for each pixel will depend on the neighborhood of the current pixel.By assuming that the pixels in the foreground and background both obey Gaussian distribution, we can construct the following expression forPIC:(7)PIC=∏m∈M12πσCmexp−Im−μCm2σCm22,where μCm and σCm are the mean and standard deviation of the pixel values in the foreground and background regions, respectively. By using the level set function ϕ and instantiating “C” as “F” and “B”, we can give their calculation formulas:(8)μFm=meanIminϕ≥0,σFm=stdIminϕ≥0,μBm=meanIminϕ<0,σBm=stdIminϕ<0.As forPC, we assume that the label field of the current pixel is only associated with its neighborhood, which is reasonable because it is highly consistent with the Markov property. Therefore, based on the Hammersley–Clifford theorem [13], PC can be expressed as a Gibbs density of the form:(9)PC=∏m∈M1Nexp−∑k∈ℕwmEkCm,where N is a constant for normalization, k is a set, which is made up of all the possible cliques, and EkCm is the clique energy function, and it has the following form of calculation formula:(10)EkCm=−αexp−Im−IpmCm=Cpm,α21+exp−Im−Ipm−1Cm≠Cpm,where pm∈ℕwm, Im−Ipm is the absolute value of the difference between the center pixel and one of its 8 neighborhood pixels and α is a Gibbsian parameter which usually needs to be set as an appropriate constant. From the calculation formula of EkCm, we can obviously find that its value is directly related to the local intensity information of the image; therefore, the prior segmentation probability PC has strong adaptability to the local image data. This is the important reason why the segmentation model of this paper can achieve excellent segmentation performance.According to the theory of optimization, we know that the MAP problem is equivalent to the energy minimization problem; therefore, we can restatePCI to the following form:(11)EMAPCI=−logPCI≈−logPC⋅PIC=−logPC−logPIC,=∑m∈MIm−μCm22σCm2+log2πZσCm+∑k∈ℕwmEkCm.By instantiating “C” as “F” and “B”, at the same time, introducing the level set function ϕ and Heaviside function H, we can construct the MRF energy functional EMRF, which has the following linear combination form:(12)EMRF=EMAPFI⋅Hϕ+EMAPBI⋅1−Hϕ.


### 3.2. Gradient Descent Flow and Graph Cuts Optimization-Based Numerical Implementation

#### 3.2.1. Gradient Descent Flow

By combining the aforementioned subenergy terms together, we obtain the total energy functional corresponding to the proposed segmentation model:(13)Eϕ=EREGϕ+EMRFϕ,=∬Ωδϕ∇ϕf∇Gσ∗If∇Gσ∗Idxdy+EMAPFI⋅Hϕ+EMAPBI⋅1−Hϕ.In order to derive the Euler–Lagrange equation corresponding to the level set functionϕ, here, we use the slightly regularized forms of the two functions δϕ and Hϕ, their formulas are as follows:(14)Hεz=121+2πarctanzε,δεz=1π⋅εε2+z2.Next, by introducing an artificial time variablet and minimizing the energy functional Eϕ with respect to ϕ, we obtain the following gradient descent flow equation:(15)∂ϕ∂t=−∂Eϕ∂ϕ=δεϕdiv∇ϕf∇Gσ∗I−2∇ϕ−1f∇Gσ∗I−1δε′∇ϕf∇Gσ∗I−δεϕEMAPFI−EMAPBI,where δε′ is the partial derivative of δεϕ with respect to ϕ. Based on this equation, we can realize the iterative updating of the level set function ϕ, and then realize the dynamic evolution of the curve.If the traditional numerical discretization methods, such as the finite difference type methods used in our previous conference paper [14], are used to solve the gradient descent equation shown in equation (15), the time consumption of the evolution process will be very large, which is not feasible for real-time class applications. In the next section, by transforming the problem here into a graph cuts optimization problem, we can design a fast and efficient numerical solution.


#### 3.2.2. Graph Cuts Optimization-Based Numerical Implementation

In this section, we will efficiently solve the proposed level set model from the perspective of graph cuts optimization. Because the graph cuts algorithms have two obvious features: rapidity (the time consumption of the iteration process is small) and globality (the iteration process does not fall into the local optimal point), the numerical implementation schemes based on such algorithm frameworks can greatly improve the speed of curve evolution process and the accuracy of segmentation results, and the final segmentation results are weakly correlated with the initial position of active contours (the benefit of the globality).SupposeG is a graph, which has an overall architecture like G=V,E,W, where V is the set of all vertices, E is the set of all connected edges, and W is the set of weights (greater than or equal to 0) associated with each connected edge. In order to intuitively express the segmentation and classification process of the graph, it is usually necessary to set two special terminal nodes: s (source node) and t (sink node). For each point gp in the vertex set G, if there is a path starting from s via gp to t, then this type of graph is called a flow network. In this type of network, those edges starting from s and those entering t are called t-link, and the edges connecting vertices other than s and t are called n-link. A cut of a graph can divide it into two subsets S and T without intersections, and we have s∈S and t∈T. For a given cut S,T, its cost function can be defined as follows:(16)CG=∑s∈S,t∈Tωst,where ωst is the weight of the connection edge of nodes s and t.How to effectively estimate the length of the contour within the graph cuts framework is a critical step for our algorithm. Fortunately, this issue has been thoroughly studied by Kolmogrov and Boykov [14–16]. They construct a grid graph and assign the weight values to the connection edges of the graph so that we can associate the cost of the cut and the length of the contour intuitively with the help of the Cauchy–Crofton formula. The Cauchy–Crofton formula explains that, by drawing sufficiently large number of straight lines in all direction from 0 to 2π and counting the number of points of intersections of the lines and the contour of interest, an estimate of the length of the contour can be obtained. Then, by some reasonably partitioning, the following discrete formula can be used to approximate the length of the contours:(17)CE=∑m=1Mncmδ2⋅Δθm2⋅em,where M is the total number of directions of the neighborhood system used. For the 4-neighborhood system shown in Figure 3, M = 2. For the 8-neighborhood and 16-neighborhood systems shown in Figure 3, their values are 4 and 8, respectively. Δθm is the angle between two adjacent directions, δ is the spacing of the digital discrete grid, em is the length of the edge em, and ncm is the total number of intersections of the active contour with all lines in the mth direction, i.e., em in Figure 3. In order to estimate the value of ncm, we introduce a metric function MFi,j to determine whether the line connecting the two pixels i and j intersects the active contour or not. It is clear that the line intersects the contour only when xi and xj have different label values. In view of this, the metric function MFi,j can be defined as follows:(18)MFi,j=xi1−xj+xj1−xi.Figure 3
Different neighborhood systems used to estimate the length of active contour. (a)ℕ4 system; (b) ℕ8 system; (c) ℕ16 system.Then, we have(19)ncm=∑i,j∈emMFi,j.If we choose the same edge weight for the straight line group with the same direction, then(20)ωm=δ2⋅Δθmem.Combining equation (17) with equation (20), the contour length CE can be rewritten as follows:(21)CE=∑m=1Mncm⋅ωm≈CG.(1) Discrete Representation of Energy Functional. The energy functional equation shown in equation (13) is defined in continuous data field. In order to associate it with the aforementioned grid graph, we first need to discretize it. Suppose GP=r,cr=1,2,…,row;c=1,2,…,col is the set of pixel grid points of the image to be considered, and the horizontal and vertical spacing between grid points is 1. Obviously, the total number of grid points in set GP is NGP=row×col. To facilitate the representation of a partition of the grid graph, we define a binary flag function xrcxrc∈0,1,r,c∈GP similar to that used in [17]. Its definition with the form of a simple expression i=r,c∈GP is as follows:(22)xi=1,ifϕi≥0,0,otherwise.Based on the abovementioned flag function, the first term of equation (13) can be discretized as the following expression:(23)E1=∑i∈GPαixi+βi1−xi,where α=EMAPFI and β=EMAPBI.Based on the same processing logic, we can obtain the discrete forms of the second to fourth terms based on the above flag function:(24)E2=∑i∈GP∑j∈ℕrjϕi−ϕjηi−1ηi⋅ωij1−xixj+1−xjxi,where η=f∇Gσ∗I.By combining the above various subitems, a total discretization form corresponding to equation (13) can be obtained in the following form:(25)Ex,ϕ=∑i∈GPαixi+βi1−xi+∑i∈GP∑j∈ℕrjϕi−ϕjηi−1ηi⋅ωij1−xixj+1−xjxi.(2) Energy Minimization Using Graph Cuts Optimization. According to the function classification rules provided in [18], it is easy to see that equation (25) is a typical function that can be optimized using the graph cuts framework. Suppose the framework of the flow network to be considered is G=V,E,W, which has a source node named s and a sink node named t. By treating the image pixels as the grid graph nodes, we can closely associate the image with the flow network G. For each pixel i, there are two t-links, which are s,i and i,t, and their weights are defined as wsi and wit, respectively. In addition, each pair pixel i,j in the neighborhood system is connected by n-link, where we assume that its connection weight is wij. By analyzing equation (25) and combining the similar items in it, we can get the formula for calculating the various weights of G as follows:(26)wsi=αi,i∈GP,(27)wit=βi,i∈GP,(28)wij=ϕi−ϕjηi−1ηiωij,i∈GP,∀j∈ℕrj.After obtaining the weight coefficients shown in equations (26)–(28), we can obtain the minimum cut of G by some typical maximum flow minimum cut algorithms [19–21]. The principles of graph cuts tell us that the minimum cut of the graph corresponds exactly to the minimum value of the energy functional shown in equation (25). After completing each iteration, we can get the updated binary label xixi∈0,1,i=1,2,…,NGP and level set function ϕ.


### 3.3. Algorithm Steps

The detailed execution steps of the proposed algorithm are as follows:(1)
Place (manually set or automatically generated via a computer program) an initial curveC on the image data to be processed and calculate the corresponding signed distance function, which is exactly the initial level set function ϕ required for the evolution process. The symbol distribution rule of ϕ is as follows: ϕx,y>0 if point x,y is inside C, ϕx,y<0 if x,y is outside C, and ϕx,y=0 if x,y is on C.(2)
Computeη according to equation (4).(3)
Computeα and β according to equation (11).(4)
Let the values ofxi,i=1,2,…,NGP corresponding to ϕx,y>0 and ϕx,y<0 be 1 and 0, respectively.(5)
Use equations (26)–(28) to construct the grid graph G and use certain graph cuts algorithms to derive the minimum cut of G until a set of updated binary labels xixi∈0,1,i=1,2,…,NGP are obtained.(6)
Repeat step (5) until the graph cuts optimization process reaches its convergence state, wherein the portion of the binary label value equal to 1 is the final segmentation result.


## 3.1. Energy Terms

### 3.1.1. Regularization Term

In order to enhance the smoothness of the level set evolution process, we introduce an adaptive regularization constraint that is directly applied to the zero level set, the effect of which is to reduce the influence of noise and ensure that the active contour does not pass through the weak object boundaries.In order to solve the problem of regularization of zero level curve, Zhou and Mu [10] proposed an efficient regularization scheme based on weighted Dirichlet integral. Although this method has achieved good results, its constant exponent gives it the following problems: cannot effectively reflect the local characteristics of the image; therefore, the regularization scheme cannot let the exponent automatically adapt to the image data. In order to enhance the adaptability of the constraint to the local image information, we propose the following weighted Dirichlet integral regularization with variable exponent:(3)EREGϕ=∬Ωδϕ∇ϕf∇Gσ∗If∇Gσ∗Idxdy,where “Gσ∗I” is the convolution operation between the image I and the Gaussian filter kernel Gσ (whose standard deviation is σ), “∇” is a gradient operator, and the exponent fn is a strictly monotone function and has the following limiting properties: (i) limn⟶0fn=1; (ii) limn⟶+∞fn=1.5. By matching these two features, we can construct the following simple function:(4)fn=1.5−1n+2.Below, we briefly analyze the macroscopic properties of the functionEREG−ZEROϕ:(i)
Whenf∇Gσ∗I equals to 1, RLSCϕ will degenerate to the following simple form: RREGϕ=∬Ωδϕ∇ϕdxdy. By comparing it with the C-V model, we find that this degenerate form is the curve length constraint used in the C-V model, which calculates the length of the zero curve of the level set function ϕ.(ii)
Whenf∇Gσ∗I is a constant F greater than 1, EREGϕ will be simplified as a geometric constraint for zero level curve used in the literature [10].For most images, the intensity value of the same foreground area is not strictly uniform, that is to say, different local areas within the foreground range have different local characteristics. Under this generality, if the above exponentf⋅ is set as a constant, the local characteristics of the image cannot be correctly reflected; in other words, the exponent f⋅ cannot automatically match the image data. Therefore, it is very unreasonable to set the exponent f⋅ constant. In view of this, we take a different approach. In our framework, the exponent f⋅ is directly related to the local intensity information of the image. Our approach ensures that weak regularization (f∇Gσ∗I≈1) is performed within regions with almost constant intensity values (the gradient value of the image is almost zero); as a result, it can effectively avoid the problem of the disappearance of weak boundaries. In contrast, in other areas, it will take a strong regularization operation to force the removal of false contours. Therefore, the existence of variable exponent f⋅ makes our regularization scheme adaptively choose between weak regularization and strong regularization. In the experimental results section, we will through a simple example to intuitively feel the effect of the variable exponent strategy.


### 3.1.2. MRF-Based External Energy Term

The traditional level set segmentation models ignore the coupling relationship between neighboring pixels, so the evolution curve often converges to the wrong noise position. In order to improve the overall antinoise performance of the level set segmentation model, we introduce the idea of MRF to model the coupling relationship between the pixels and their neighborhoods.Firstly, we represent the square neighborhood required by the MRF idea asℕwm, whose size is 2w+1×2w+1 and m refers to the component of the random variable M. In order to match the requirements of MRF, we need to first construct a coupled random field whose purpose is to design our variable and symbolic system, here we represent the coupled random field as CRF=I,C, where I=Im,m∈M represents the field of pixel intensity, and C=Cm,m∈M represents the label random field (also known as characteristic random field), which is used to distinguish the foreground (F) from the background (B).In addition, the relevant theoretical knowledge of MFR tells us that the energy functional based on MRF finds the optimal category (foreground or background) label by using the neighborhood information of each pixel. Usually, we can achieve this goal by maximizing a posteriori (MAP) probability distribution function. With the help of Bayes theory, we can expressPCI as follows:(5)PCI=PIC⋅PCPI,where PCI, PC, PIC, and PI represent segmentation probability, priori segmentation probability, conditional segmentation probability, and image probability, respectively. When the image data are given, the value of PI is a constant; therefore, we can delete it from the above expression, and the deletion will not affect the proportional relationship contained in the equation. Therefore, we have the following further proportional relationship:(6)PCI∝PIC⋅PC.The existing level set segmentation models generally assume that the pixels in the foreground and background regions are independent from each other. However, such assumption is not reasonable and directly leads to their segmentation results vulnerable to noise interference. In order to overcome this shortcoming, we introduce the MRF theory to model the two components ofPIC and PC in the above expression, that is to say, when we construct the probability expression, we incorporate the neighborhood relationships into our consideration system. After the introduction of MRF, the values of PIC and PC for each pixel will depend on the neighborhood of the current pixel.By assuming that the pixels in the foreground and background both obey Gaussian distribution, we can construct the following expression forPIC:(7)PIC=∏m∈M12πσCmexp−Im−μCm2σCm22,where μCm and σCm are the mean and standard deviation of the pixel values in the foreground and background regions, respectively. By using the level set function ϕ and instantiating “C” as “F” and “B”, we can give their calculation formulas:(8)μFm=meanIminϕ≥0,σFm=stdIminϕ≥0,μBm=meanIminϕ<0,σBm=stdIminϕ<0.As forPC, we assume that the label field of the current pixel is only associated with its neighborhood, which is reasonable because it is highly consistent with the Markov property. Therefore, based on the Hammersley–Clifford theorem [13], PC can be expressed as a Gibbs density of the form:(9)PC=∏m∈M1Nexp−∑k∈ℕwmEkCm,where N is a constant for normalization, k is a set, which is made up of all the possible cliques, and EkCm is the clique energy function, and it has the following form of calculation formula:(10)EkCm=−αexp−Im−IpmCm=Cpm,α21+exp−Im−Ipm−1Cm≠Cpm,where pm∈ℕwm, Im−Ipm is the absolute value of the difference between the center pixel and one of its 8 neighborhood pixels and α is a Gibbsian parameter which usually needs to be set as an appropriate constant. From the calculation formula of EkCm, we can obviously find that its value is directly related to the local intensity information of the image; therefore, the prior segmentation probability PC has strong adaptability to the local image data. This is the important reason why the segmentation model of this paper can achieve excellent segmentation performance.According to the theory of optimization, we know that the MAP problem is equivalent to the energy minimization problem; therefore, we can restatePCI to the following form:(11)EMAPCI=−logPCI≈−logPC⋅PIC=−logPC−logPIC,=∑m∈MIm−μCm22σCm2+log2πZσCm+∑k∈ℕwmEkCm.By instantiating “C” as “F” and “B”, at the same time, introducing the level set function ϕ and Heaviside function H, we can construct the MRF energy functional EMRF, which has the following linear combination form:(12)EMRF=EMAPFI⋅Hϕ+EMAPBI⋅1−Hϕ.


## 3.1.1. Regularization Term

In order to enhance the smoothness of the level set evolution process, we introduce an adaptive regularization constraint that is directly applied to the zero level set, the effect of which is to reduce the influence of noise and ensure that the active contour does not pass through the weak object boundaries.In order to solve the problem of regularization of zero level curve, Zhou and Mu [10] proposed an efficient regularization scheme based on weighted Dirichlet integral. Although this method has achieved good results, its constant exponent gives it the following problems: cannot effectively reflect the local characteristics of the image; therefore, the regularization scheme cannot let the exponent automatically adapt to the image data. In order to enhance the adaptability of the constraint to the local image information, we propose the following weighted Dirichlet integral regularization with variable exponent:(3)EREGϕ=∬Ωδϕ∇ϕf∇Gσ∗If∇Gσ∗Idxdy,where “Gσ∗I” is the convolution operation between the image I and the Gaussian filter kernel Gσ (whose standard deviation is σ), “∇” is a gradient operator, and the exponent fn is a strictly monotone function and has the following limiting properties: (i) limn⟶0fn=1; (ii) limn⟶+∞fn=1.5. By matching these two features, we can construct the following simple function:(4)fn=1.5−1n+2.Below, we briefly analyze the macroscopic properties of the functionEREG−ZEROϕ:(i)
Whenf∇Gσ∗I equals to 1, RLSCϕ will degenerate to the following simple form: RREGϕ=∬Ωδϕ∇ϕdxdy. By comparing it with the C-V model, we find that this degenerate form is the curve length constraint used in the C-V model, which calculates the length of the zero curve of the level set function ϕ.(ii)
Whenf∇Gσ∗I is a constant F greater than 1, EREGϕ will be simplified as a geometric constraint for zero level curve used in the literature [10].For most images, the intensity value of the same foreground area is not strictly uniform, that is to say, different local areas within the foreground range have different local characteristics. Under this generality, if the above exponentf⋅ is set as a constant, the local characteristics of the image cannot be correctly reflected; in other words, the exponent f⋅ cannot automatically match the image data. Therefore, it is very unreasonable to set the exponent f⋅ constant. In view of this, we take a different approach. In our framework, the exponent f⋅ is directly related to the local intensity information of the image. Our approach ensures that weak regularization (f∇Gσ∗I≈1) is performed within regions with almost constant intensity values (the gradient value of the image is almost zero); as a result, it can effectively avoid the problem of the disappearance of weak boundaries. In contrast, in other areas, it will take a strong regularization operation to force the removal of false contours. Therefore, the existence of variable exponent f⋅ makes our regularization scheme adaptively choose between weak regularization and strong regularization. In the experimental results section, we will through a simple example to intuitively feel the effect of the variable exponent strategy.


## 3.1.2. MRF-Based External Energy Term

The traditional level set segmentation models ignore the coupling relationship between neighboring pixels, so the evolution curve often converges to the wrong noise position. In order to improve the overall antinoise performance of the level set segmentation model, we introduce the idea of MRF to model the coupling relationship between the pixels and their neighborhoods.Firstly, we represent the square neighborhood required by the MRF idea asℕwm, whose size is 2w+1×2w+1 and m refers to the component of the random variable M. In order to match the requirements of MRF, we need to first construct a coupled random field whose purpose is to design our variable and symbolic system, here we represent the coupled random field as CRF=I,C, where I=Im,m∈M represents the field of pixel intensity, and C=Cm,m∈M represents the label random field (also known as characteristic random field), which is used to distinguish the foreground (F) from the background (B).In addition, the relevant theoretical knowledge of MFR tells us that the energy functional based on MRF finds the optimal category (foreground or background) label by using the neighborhood information of each pixel. Usually, we can achieve this goal by maximizing a posteriori (MAP) probability distribution function. With the help of Bayes theory, we can expressPCI as follows:(5)PCI=PIC⋅PCPI,where PCI, PC, PIC, and PI represent segmentation probability, priori segmentation probability, conditional segmentation probability, and image probability, respectively. When the image data are given, the value of PI is a constant; therefore, we can delete it from the above expression, and the deletion will not affect the proportional relationship contained in the equation. Therefore, we have the following further proportional relationship:(6)PCI∝PIC⋅PC.The existing level set segmentation models generally assume that the pixels in the foreground and background regions are independent from each other. However, such assumption is not reasonable and directly leads to their segmentation results vulnerable to noise interference. In order to overcome this shortcoming, we introduce the MRF theory to model the two components ofPIC and PC in the above expression, that is to say, when we construct the probability expression, we incorporate the neighborhood relationships into our consideration system. After the introduction of MRF, the values of PIC and PC for each pixel will depend on the neighborhood of the current pixel.By assuming that the pixels in the foreground and background both obey Gaussian distribution, we can construct the following expression forPIC:(7)PIC=∏m∈M12πσCmexp−Im−μCm2σCm22,where μCm and σCm are the mean and standard deviation of the pixel values in the foreground and background regions, respectively. By using the level set function ϕ and instantiating “C” as “F” and “B”, we can give their calculation formulas:(8)μFm=meanIminϕ≥0,σFm=stdIminϕ≥0,μBm=meanIminϕ<0,σBm=stdIminϕ<0.As forPC, we assume that the label field of the current pixel is only associated with its neighborhood, which is reasonable because it is highly consistent with the Markov property. Therefore, based on the Hammersley–Clifford theorem [13], PC can be expressed as a Gibbs density of the form:(9)PC=∏m∈M1Nexp−∑k∈ℕwmEkCm,where N is a constant for normalization, k is a set, which is made up of all the possible cliques, and EkCm is the clique energy function, and it has the following form of calculation formula:(10)EkCm=−αexp−Im−IpmCm=Cpm,α21+exp−Im−Ipm−1Cm≠Cpm,where pm∈ℕwm, Im−Ipm is the absolute value of the difference between the center pixel and one of its 8 neighborhood pixels and α is a Gibbsian parameter which usually needs to be set as an appropriate constant. From the calculation formula of EkCm, we can obviously find that its value is directly related to the local intensity information of the image; therefore, the prior segmentation probability PC has strong adaptability to the local image data. This is the important reason why the segmentation model of this paper can achieve excellent segmentation performance.According to the theory of optimization, we know that the MAP problem is equivalent to the energy minimization problem; therefore, we can restatePCI to the following form:(11)EMAPCI=−logPCI≈−logPC⋅PIC=−logPC−logPIC,=∑m∈MIm−μCm22σCm2+log2πZσCm+∑k∈ℕwmEkCm.By instantiating “C” as “F” and “B”, at the same time, introducing the level set function ϕ and Heaviside function H, we can construct the MRF energy functional EMRF, which has the following linear combination form:(12)EMRF=EMAPFI⋅Hϕ+EMAPBI⋅1−Hϕ.


## 3.2. Gradient Descent Flow and Graph Cuts Optimization-Based Numerical Implementation

### 3.2.1. Gradient Descent Flow

By combining the aforementioned subenergy terms together, we obtain the total energy functional corresponding to the proposed segmentation model:(13)Eϕ=EREGϕ+EMRFϕ,=∬Ωδϕ∇ϕf∇Gσ∗If∇Gσ∗Idxdy+EMAPFI⋅Hϕ+EMAPBI⋅1−Hϕ.In order to derive the Euler–Lagrange equation corresponding to the level set functionϕ, here, we use the slightly regularized forms of the two functions δϕ and Hϕ, their formulas are as follows:(14)Hεz=121+2πarctanzε,δεz=1π⋅εε2+z2.Next, by introducing an artificial time variablet and minimizing the energy functional Eϕ with respect to ϕ, we obtain the following gradient descent flow equation:(15)∂ϕ∂t=−∂Eϕ∂ϕ=δεϕdiv∇ϕf∇Gσ∗I−2∇ϕ−1f∇Gσ∗I−1δε′∇ϕf∇Gσ∗I−δεϕEMAPFI−EMAPBI,where δε′ is the partial derivative of δεϕ with respect to ϕ. Based on this equation, we can realize the iterative updating of the level set function ϕ, and then realize the dynamic evolution of the curve.If the traditional numerical discretization methods, such as the finite difference type methods used in our previous conference paper [14], are used to solve the gradient descent equation shown in equation (15), the time consumption of the evolution process will be very large, which is not feasible for real-time class applications. In the next section, by transforming the problem here into a graph cuts optimization problem, we can design a fast and efficient numerical solution.


### 3.2.2. Graph Cuts Optimization-Based Numerical Implementation

In this section, we will efficiently solve the proposed level set model from the perspective of graph cuts optimization. Because the graph cuts algorithms have two obvious features: rapidity (the time consumption of the iteration process is small) and globality (the iteration process does not fall into the local optimal point), the numerical implementation schemes based on such algorithm frameworks can greatly improve the speed of curve evolution process and the accuracy of segmentation results, and the final segmentation results are weakly correlated with the initial position of active contours (the benefit of the globality).SupposeG is a graph, which has an overall architecture like G=V,E,W, where V is the set of all vertices, E is the set of all connected edges, and W is the set of weights (greater than or equal to 0) associated with each connected edge. In order to intuitively express the segmentation and classification process of the graph, it is usually necessary to set two special terminal nodes: s (source node) and t (sink node). For each point gp in the vertex set G, if there is a path starting from s via gp to t, then this type of graph is called a flow network. In this type of network, those edges starting from s and those entering t are called t-link, and the edges connecting vertices other than s and t are called n-link. A cut of a graph can divide it into two subsets S and T without intersections, and we have s∈S and t∈T. For a given cut S,T, its cost function can be defined as follows:(16)CG=∑s∈S,t∈Tωst,where ωst is the weight of the connection edge of nodes s and t.How to effectively estimate the length of the contour within the graph cuts framework is a critical step for our algorithm. Fortunately, this issue has been thoroughly studied by Kolmogrov and Boykov [14–16]. They construct a grid graph and assign the weight values to the connection edges of the graph so that we can associate the cost of the cut and the length of the contour intuitively with the help of the Cauchy–Crofton formula. The Cauchy–Crofton formula explains that, by drawing sufficiently large number of straight lines in all direction from 0 to 2π and counting the number of points of intersections of the lines and the contour of interest, an estimate of the length of the contour can be obtained. Then, by some reasonably partitioning, the following discrete formula can be used to approximate the length of the contours:(17)CE=∑m=1Mncmδ2⋅Δθm2⋅em,where M is the total number of directions of the neighborhood system used. For the 4-neighborhood system shown in Figure 3, M = 2. For the 8-neighborhood and 16-neighborhood systems shown in Figure 3, their values are 4 and 8, respectively. Δθm is the angle between two adjacent directions, δ is the spacing of the digital discrete grid, em is the length of the edge em, and ncm is the total number of intersections of the active contour with all lines in the mth direction, i.e., em in Figure 3. In order to estimate the value of ncm, we introduce a metric function MFi,j to determine whether the line connecting the two pixels i and j intersects the active contour or not. It is clear that the line intersects the contour only when xi and xj have different label values. In view of this, the metric function MFi,j can be defined as follows:(18)MFi,j=xi1−xj+xj1−xi.Figure 3
Different neighborhood systems used to estimate the length of active contour. (a)ℕ4 system; (b) ℕ8 system; (c) ℕ16 system.Then, we have(19)ncm=∑i,j∈emMFi,j.If we choose the same edge weight for the straight line group with the same direction, then(20)ωm=δ2⋅Δθmem.Combining equation (17) with equation (20), the contour length CE can be rewritten as follows:(21)CE=∑m=1Mncm⋅ωm≈CG.(1) Discrete Representation of Energy Functional. The energy functional equation shown in equation (13) is defined in continuous data field. In order to associate it with the aforementioned grid graph, we first need to discretize it. Suppose GP=r,cr=1,2,…,row;c=1,2,…,col is the set of pixel grid points of the image to be considered, and the horizontal and vertical spacing between grid points is 1. Obviously, the total number of grid points in set GP is NGP=row×col. To facilitate the representation of a partition of the grid graph, we define a binary flag function xrcxrc∈0,1,r,c∈GP similar to that used in [17]. Its definition with the form of a simple expression i=r,c∈GP is as follows:(22)xi=1,ifϕi≥0,0,otherwise.Based on the abovementioned flag function, the first term of equation (13) can be discretized as the following expression:(23)E1=∑i∈GPαixi+βi1−xi,where α=EMAPFI and β=EMAPBI.Based on the same processing logic, we can obtain the discrete forms of the second to fourth terms based on the above flag function:(24)E2=∑i∈GP∑j∈ℕrjϕi−ϕjηi−1ηi⋅ωij1−xixj+1−xjxi,where η=f∇Gσ∗I.By combining the above various subitems, a total discretization form corresponding to equation (13) can be obtained in the following form:(25)Ex,ϕ=∑i∈GPαixi+βi1−xi+∑i∈GP∑j∈ℕrjϕi−ϕjηi−1ηi⋅ωij1−xixj+1−xjxi.(2) Energy Minimization Using Graph Cuts Optimization. According to the function classification rules provided in [18], it is easy to see that equation (25) is a typical function that can be optimized using the graph cuts framework. Suppose the framework of the flow network to be considered is G=V,E,W, which has a source node named s and a sink node named t. By treating the image pixels as the grid graph nodes, we can closely associate the image with the flow network G. For each pixel i, there are two t-links, which are s,i and i,t, and their weights are defined as wsi and wit, respectively. In addition, each pair pixel i,j in the neighborhood system is connected by n-link, where we assume that its connection weight is wij. By analyzing equation (25) and combining the similar items in it, we can get the formula for calculating the various weights of G as follows:(26)wsi=αi,i∈GP,(27)wit=βi,i∈GP,(28)wij=ϕi−ϕjηi−1ηiωij,i∈GP,∀j∈ℕrj.After obtaining the weight coefficients shown in equations (26)–(28), we can obtain the minimum cut of G by some typical maximum flow minimum cut algorithms [19–21]. The principles of graph cuts tell us that the minimum cut of the graph corresponds exactly to the minimum value of the energy functional shown in equation (25). After completing each iteration, we can get the updated binary label xixi∈0,1,i=1,2,…,NGP and level set function ϕ.


## 3.2.1. Gradient Descent Flow

By combining the aforementioned subenergy terms together, we obtain the total energy functional corresponding to the proposed segmentation model:(13)Eϕ=EREGϕ+EMRFϕ,=∬Ωδϕ∇ϕf∇Gσ∗If∇Gσ∗Idxdy+EMAPFI⋅Hϕ+EMAPBI⋅1−Hϕ.In order to derive the Euler–Lagrange equation corresponding to the level set functionϕ, here, we use the slightly regularized forms of the two functions δϕ and Hϕ, their formulas are as follows:(14)Hεz=121+2πarctanzε,δεz=1π⋅εε2+z2.Next, by introducing an artificial time variablet and minimizing the energy functional Eϕ with respect to ϕ, we obtain the following gradient descent flow equation:(15)∂ϕ∂t=−∂Eϕ∂ϕ=δεϕdiv∇ϕf∇Gσ∗I−2∇ϕ−1f∇Gσ∗I−1δε′∇ϕf∇Gσ∗I−δεϕEMAPFI−EMAPBI,where δε′ is the partial derivative of δεϕ with respect to ϕ. Based on this equation, we can realize the iterative updating of the level set function ϕ, and then realize the dynamic evolution of the curve.If the traditional numerical discretization methods, such as the finite difference type methods used in our previous conference paper [14], are used to solve the gradient descent equation shown in equation (15), the time consumption of the evolution process will be very large, which is not feasible for real-time class applications. In the next section, by transforming the problem here into a graph cuts optimization problem, we can design a fast and efficient numerical solution.


## 3.2.2. Graph Cuts Optimization-Based Numerical Implementation

In this section, we will efficiently solve the proposed level set model from the perspective of graph cuts optimization. Because the graph cuts algorithms have two obvious features: rapidity (the time consumption of the iteration process is small) and globality (the iteration process does not fall into the local optimal point), the numerical implementation schemes based on such algorithm frameworks can greatly improve the speed of curve evolution process and the accuracy of segmentation results, and the final segmentation results are weakly correlated with the initial position of active contours (the benefit of the globality).SupposeG is a graph, which has an overall architecture like G=V,E,W, where V is the set of all vertices, E is the set of all connected edges, and W is the set of weights (greater than or equal to 0) associated with each connected edge. In order to intuitively express the segmentation and classification process of the graph, it is usually necessary to set two special terminal nodes: s (source node) and t (sink node). For each point gp in the vertex set G, if there is a path starting from s via gp to t, then this type of graph is called a flow network. In this type of network, those edges starting from s and those entering t are called t-link, and the edges connecting vertices other than s and t are called n-link. A cut of a graph can divide it into two subsets S and T without intersections, and we have s∈S and t∈T. For a given cut S,T, its cost function can be defined as follows:(16)CG=∑s∈S,t∈Tωst,where ωst is the weight of the connection edge of nodes s and t.How to effectively estimate the length of the contour within the graph cuts framework is a critical step for our algorithm. Fortunately, this issue has been thoroughly studied by Kolmogrov and Boykov [14–16]. They construct a grid graph and assign the weight values to the connection edges of the graph so that we can associate the cost of the cut and the length of the contour intuitively with the help of the Cauchy–Crofton formula. The Cauchy–Crofton formula explains that, by drawing sufficiently large number of straight lines in all direction from 0 to 2π and counting the number of points of intersections of the lines and the contour of interest, an estimate of the length of the contour can be obtained. Then, by some reasonably partitioning, the following discrete formula can be used to approximate the length of the contours:(17)CE=∑m=1Mncmδ2⋅Δθm2⋅em,where M is the total number of directions of the neighborhood system used. For the 4-neighborhood system shown in Figure 3, M = 2. For the 8-neighborhood and 16-neighborhood systems shown in Figure 3, their values are 4 and 8, respectively. Δθm is the angle between two adjacent directions, δ is the spacing of the digital discrete grid, em is the length of the edge em, and ncm is the total number of intersections of the active contour with all lines in the mth direction, i.e., em in Figure 3. In order to estimate the value of ncm, we introduce a metric function MFi,j to determine whether the line connecting the two pixels i and j intersects the active contour or not. It is clear that the line intersects the contour only when xi and xj have different label values. In view of this, the metric function MFi,j can be defined as follows:(18)MFi,j=xi1−xj+xj1−xi.Figure 3
Different neighborhood systems used to estimate the length of active contour. (a)ℕ4 system; (b) ℕ8 system; (c) ℕ16 system.Then, we have(19)ncm=∑i,j∈emMFi,j.If we choose the same edge weight for the straight line group with the same direction, then(20)ωm=δ2⋅Δθmem.Combining equation (17) with equation (20), the contour length CE can be rewritten as follows:(21)CE=∑m=1Mncm⋅ωm≈CG.(1) Discrete Representation of Energy Functional. The energy functional equation shown in equation (13) is defined in continuous data field. In order to associate it with the aforementioned grid graph, we first need to discretize it. Suppose GP=r,cr=1,2,…,row;c=1,2,…,col is the set of pixel grid points of the image to be considered, and the horizontal and vertical spacing between grid points is 1. Obviously, the total number of grid points in set GP is NGP=row×col. To facilitate the representation of a partition of the grid graph, we define a binary flag function xrcxrc∈0,1,r,c∈GP similar to that used in [17]. Its definition with the form of a simple expression i=r,c∈GP is as follows:(22)xi=1,ifϕi≥0,0,otherwise.Based on the abovementioned flag function, the first term of equation (13) can be discretized as the following expression:(23)E1=∑i∈GPαixi+βi1−xi,where α=EMAPFI and β=EMAPBI.Based on the same processing logic, we can obtain the discrete forms of the second to fourth terms based on the above flag function:(24)E2=∑i∈GP∑j∈ℕrjϕi−ϕjηi−1ηi⋅ωij1−xixj+1−xjxi,where η=f∇Gσ∗I.By combining the above various subitems, a total discretization form corresponding to equation (13) can be obtained in the following form:(25)Ex,ϕ=∑i∈GPαixi+βi1−xi+∑i∈GP∑j∈ℕrjϕi−ϕjηi−1ηi⋅ωij1−xixj+1−xjxi.(2) Energy Minimization Using Graph Cuts Optimization. According to the function classification rules provided in [18], it is easy to see that equation (25) is a typical function that can be optimized using the graph cuts framework. Suppose the framework of the flow network to be considered is G=V,E,W, which has a source node named s and a sink node named t. By treating the image pixels as the grid graph nodes, we can closely associate the image with the flow network G. For each pixel i, there are two t-links, which are s,i and i,t, and their weights are defined as wsi and wit, respectively. In addition, each pair pixel i,j in the neighborhood system is connected by n-link, where we assume that its connection weight is wij. By analyzing equation (25) and combining the similar items in it, we can get the formula for calculating the various weights of G as follows:(26)wsi=αi,i∈GP,(27)wit=βi,i∈GP,(28)wij=ϕi−ϕjηi−1ηiωij,i∈GP,∀j∈ℕrj.After obtaining the weight coefficients shown in equations (26)–(28), we can obtain the minimum cut of G by some typical maximum flow minimum cut algorithms [19–21]. The principles of graph cuts tell us that the minimum cut of the graph corresponds exactly to the minimum value of the energy functional shown in equation (25). After completing each iteration, we can get the updated binary label xixi∈0,1,i=1,2,…,NGP and level set function ϕ.


## 3.3. Algorithm Steps

The detailed execution steps of the proposed algorithm are as follows:(1)
Place (manually set or automatically generated via a computer program) an initial curveC on the image data to be processed and calculate the corresponding signed distance function, which is exactly the initial level set function ϕ required for the evolution process. The symbol distribution rule of ϕ is as follows: ϕx,y>0 if point x,y is inside C, ϕx,y<0 if x,y is outside C, and ϕx,y=0 if x,y is on C.(2)
Computeη according to equation (4).(3)
Computeα and β according to equation (11).(4)
Let the values ofxi,i=1,2,…,NGP corresponding to ϕx,y>0 and ϕx,y<0 be 1 and 0, respectively.(5)
Use equations (26)–(28) to construct the grid graph G and use certain graph cuts algorithms to derive the minimum cut of G until a set of updated binary labels xixi∈0,1,i=1,2,…,NGP are obtained.(6)
Repeat step (5) until the graph cuts optimization process reaches its convergence state, wherein the portion of the binary label value equal to 1 is the final segmentation result.


## 4. Experimental Results and Discussions

In this section, we shall demonstrate the experimental results of our MAR and adaptive regularization constraint-embedded level set algorithm for image segmentation. The experiments are implemented by Matlab R2012a on a computer with 2.3G Intel Core i7 CPU, 8G RAM, and Windows 7 operating system. For the following parameters, we use the same values for all experiments, i.e.,Δt=2.5, h=1, and ε=1. The standard deviation σ of the Gaussian filter kernel will change from image to image, that is to say, its value with a certain degree of image dependency. In addition, the parameter κ in the regularization term EREG−LSF follows the same change rule. In order to rigorously test and validate the proposed model, we will introduce the following methods to compare with it: GAC [1], LSR [2], CV [3], FTC [4], LBF [5], LAC [6], LPFE [7], NFT [8]; we have outlined them in Introduction. When conducting comparison experiments, we follow the same rule that all models are data-driven and do not incorporate any type of prior knowledge. In addition, all methods use the same initial curve, and the control parameters of the relevant segmentation model are set strictly according to the selection criteria given in the corresponding literature.In order to test the generalization ability of the proposed model, we adopted a method of randomly selecting the test images. Some of these images come from the literatures on image segmentation, and some come from the internet. These images are weakly related to each other.Here, we will use four region overlap metrics to compare the performances of the models quantitatively, their definitions are as follows:(29)JS=NSreference∩StestNSreference∪Stestjaccard similarity,DSC=2NSreference∩StestNSreference+NStestdice similarity coefficient,FPR=NSreference/ScommonNSreferencefalse positive ratio,FNR=NStest/ScommonNStestfalse negative ratio,where “∩” and “∪” represent the intersection and union of two regions, respectively, Stest, Sreference, and Scommon are the output region of the segmentation algorithm, the ground-truth, and the common part of two regions, respectively, N⋅ represents the number of pixels in the enclosed set. Obviously, the closer the JS [22] and DSC [23] values to 1, and the FPR and FNR values to 0, the better the segmentation results.In order to test the antinoise performance of the proposed model, we apply it to a set of noisy images, which are obtained by adding Gaussian noise of different intensities (the set of variance is0.01,0.05,0.1) to a clean synthetic image of size 128×128 (located at the upper left corner of Figure 4). Figure 4 shows the two-valued segmentation results of different methods, where the first column is the input images and the initial contours required for the level set-based segmentation methods, and the second to the ninth columns are the segmentation results by using GAC, LSR, CV, FTC, LBF, LAC, LPFE, NFT, and our method, respectively, where the red and black regions represent the foreground and background, respectively. It can be seen from the experimental results that the superiority of the proposed method in terms of visual segmentation accuracy is very obvious compared with the methods of GAC, LSR, CV, FTC, LBF and LAC, LPFE, NFT, and our method outputs the correct segmentation results at all noise levels. In the absence of noise, all methods involved in the comparisons output the correct segmentation results; in the case of noisy interference, only the LBF model outputs the correct result at a noise level of 0.01; beyond that, their outputs are either completely wrong or contain a lot of noise. Here, we will simultaneously give the corresponding quantitative assessment indicators, which are calculated by equations (24)–(27) and shown in Table 1. By comparing these quantitative data, we find that the segmentation performance of our method is optimal.Figure 4
The comparison results on a set of noisy synthetic images with increasing intensity of Gaussian noise. The first column: input image with initial contour, from the top row to the bottom row, the variances of Gaussian noise are 0, 0.01, 0.05, and 0.1, respectively. The second to the ninth columns are the segmentation results by using GAC, LSR, CV, FTC, LBF, LAC, LPFE, NFT, and our method, respectively. The parameters used in this group of experiments areσ=1.5 and κ=1.3.Table 1
Region overlap metrics of different algorithms corresponding to Figure4 numbered in order from up to down.
Image no.
Algorithms
JS
DSC
FPR
FNR
Image #1
GAC
1.0000
1.0000
0.0000
0.0000
LSR
0.9573
0.9782
0.0168
0.0268
CV
1.0000
1.0000
0.0000
0.0000
FTC
1.0000
1.0000
0.0000
0.0000
LBF
0.9946
0.9973
0.0005
0.0048
LAC
0.9800
0.9899
0.0200
0.0060
LPFE
0.9517
0.9753
0.0184
0.0310
NFT
0.9763
0.9880
0.0200
0.0039
Our method
1.0000
1.0000
0.0000
0.0000
Image #2
GAC
0.8476
0.9175
0.0671
0.0974
LSR
0.7789
0.8757
0.1499
0.0971
CV
0.8807
0.9366
0.1169
0.0031
FTC
0.8544
0.9215
0.0127
0.1436
LBF
0.9741
0.9770
0.0116
0.0043
LAC
0.9200
0.9583
0.0725
0.0087
LPFE
0.9535
0.9762
0.0244
0.0233
NFT
0.7703
0.8703
0.0308
0.2103
Our method
0.9833
0.9916
0.0114
0.0030
Image #3
GAC
0.7715
0.8710
0.1851
0.0646
LSR
0.4215
0.5931
0.1396
0.5475
CV
0.4475
0.6183
0.2917
0.4514
FTC
0.5055
0.6716
0.0119
0.4914
LBF
0.6990
0.8228
0.0087
0.2967
LAC
0.7712
0.8709
0.2045
0.0380
LPFE
0.9023
0.9487
0.0152
0.0850
NFT
0.4719
0.6412
0.0189
0.5238
Our method
0.9711
0.9853
0.0113
0.0120
Image #4
GAC
0.6547
0.7913
0.3362
0.0205
LSR
0.3446
0.5125
0.0818
0.6445
CV
0.2670
0.4214
0.3271
0.6932
FTC
0.3232
0.4886
0.0463
0.6754
LBF
0.2622
0.4154
0.0464
0.7369
LAC
0.7049
0.8269
0.2582
0.0659
LPFE
0.6453
0.7844
0.0471
0.3474
NFT
0.3689
0.5390
0.0484
0.6279
Our method
0.9429
0.9706
0.0460
0.0122Next, we apply our model to the task of medical image segmentation. It is well known that medical image segmentation is the basic task in the medical image analysis field. Accurate, robust, and rapid segmentation of the target area is an important guarantee for the subsequent quantitative analysis, three-dimensional visualization, etc.; it also lays a solid foundation for clinical applications such as image-guided surgery, radiotherapy planning, and treatment evaluation. However, medical image data generally contain noise interference, which undoubtedly poses high requirements for the robustness of segmentation algorithms. Here, we apply different level set segmentation methods to a set of medical images that come from several imaging bands. In addition to the noise interference mentioned above, there is a common feature of these images, that is, their target contours are relatively weak. The parameters associated with this set of experiments are as follows:σ=3.5 and κ=2.5. A common characteristic of the initial curves in this set of experiments is that the target areas are completely located within the large rectangular curves, so there must be a lot of background information within the active contours, and this puts a high demand on the background suppression ability of the segmentation model; in addition, for those segmentation models that are related to the initial properties (position, shape, etc.) of the curve, this type of initialization will inevitably affect their output accuracy. The comparison segmentation results are shown in Figure 5, and the corresponding quantitative assessment indicators calculated by equations (11) and (12) are shown in Table 2. From the visual comparison results shown in Figure 5, we can clearly see that the proposed model outputs perfect segmentation results on all medical images, while the segmentation results of the methods involved in the comparison are either undersegmentation or oversegmentation, or the target area is divided into a large number of fragments; in short, they all output the wrong segmentation results. Through the quantitative parameters shown in Table 2, we can further confirm the aforementioned visual comparison conclusions.Figure 5
Comparison experiments on a group of medical images with different methods. The first column: input image along with initial contours; the second column: segmentation results of GAC model; the third column: segmentation results of LSR model; the fourth column: segmentation results of CV model; the fifth column: segmentation results of FTC model; the sixth column: segmentation results of LBF model; the seventh column: segmentation results of LAC model; the eighth column: segmentation results of LPFE model; the ninth column: segmentation results of NFT model; the tenth column: segmentation results of our model.Table 2
Region overlap metrics of different algorithms corresponding to Figure5 in the same order.
Image no.
Algorithms
JS
DSC
FPR
FNR
Image #1
GAC
0.5646
0.7217
0.3699
0.1554
LSR
0.2585
0.4107
0.7415
0.3009
CV
0.7386
0.8496
0.0983
0.1967
FTC
0.7041
0.8264
0.0442
0.2722
LBF
0.6859
0.8137
0.0048
0.3119
LAC
0.8359
0.9106
0.0383
0.1353
LPFE
0.7004
0.8238
0.0418
0.2776
NFT
0.3399
0.5073
0.5828
0.3528
Our method
0.9567
0.9779
0.0035
0.0208
Image #2
GAC
0.8283
0.9061
0.1601
0.0164
LSR
0.7700
0.8700
0.1302
0.1297
CV
0.8797
0.9360
0.0162
0.1074
FTC
0.8646
0.9274
0.0202
0.1197
LBF
0.7300
0.8440
0.0200
0.2589
LAC
0.6759
0.8066
0.1711
0.2146
LPFE
0.8626
0.9262
0.0206
0.1215
NFT
0.6217
0.7667
0.2650
0.1987
Our method
0.9609
0.9800
0.0104
0.0155
Image #3
GAC
0.6232
0.7679
0.3101
0.1344
LSR
0.6243
0.7687
0.3604
0.0369
CV
0.6428
0.7826
0.2178
0.2171
FTC
0.6327
0.7750
0.2554
0.1920
LBF
0.1589
0.2743
0.3766
0.8242
LAC
0.0581
0.1099
0.7854
0.9261
LPFE
0.3726
0.5429
0.1091
0.6096
NFT
0.4171
0.5887
0.5100
0.2628
Our method
0.9494
0.9740
0.0328
0.0190
Image #4
GAC
0.7018
0.8247
0.1245
0.2204
LSR
0.5989
0.7492
0.1311
0.3416
CV
0.3637
0.5334
0.3006
0.6363
FTC
0.4062
0.5777
0.4030
0.5938
LBF
0.2854
0.4440
0.5101
0.7146
LAC
0.2692
0.4242
0.0794
0.7244
LPFE
0.4198
0.5913
0.0236
0.5800
NFT
0.5879
0.7405
0.1243
0.3586
Our method
0.9553
0.9772
0.0230
0.0227
Image #5
GAC
0.5390
0.7004
0.0532
0.4442
LSR
0.6471
0.7857
0.3051
0.0961
CV
0.5587
0.7169
0.3237
0.2374
FTC
0.5463
0.7066
0.3523
0.2226
LBF
0.5538
0.7129
0.0100
0.4430
LAC
0.5972
0.7478
0.3082
0.1864
LPFE
0.5387
0.7002
0.3597
0.2276
NFT
0.4304
0.6018
0.2255
0.5079
Our method
0.9702
0.9849
0.0081
0.0121In order to further verify the segmentation performance of the proposed model, we apply the proposed model to the task of segmentation of real infrared sea surface targets. It is well known that the infrared images at sea background are notorious for the low signal-to-noise ratio, therefore, it is often very difficult to segment infrared sea surface targets based on traditional methods. Figure6 shows a summary of the segmentation results for all the methods involved in the comparison, and Table 3 shows the corresponding quantitative indicators. When performing the initialization operation, we use the segmentation contour of OTSU [24] algorithm as the initial curve of the evolution process of level set. From this set of segmentation experiments, we can clearly see that the proposed method outputs the best segmentation results on all infrared sea surface target images; the reason is that the MRF-based level set model and the adaptive regularization strategy for the zero level curve can form a strong background suppression force; under their joint action, the proposed model can completely suppress the background interference. However, the other methods involved in the comparison do not effectively model the neighborhood pixels and the regularization strategies adopted are simple, thus they output very mediocre segmentation results.Figure 6
Comparison of different methods for infrared sea surface target image segmentation. The first column: input image; the second column: segmentation results of GAC model; the third column: segmentation results of LSR model; the fourth column: segmentation results of CV model; the fifth column: segmentation results of FTC model; the sixth column: segmentation results of LBF model; the seventh column: segmentation results of LAC model; the eighth column: segmentation results of LPFE model; the ninth column: segmentation results of NFT model; the tenth column: segmentation results of our model. We select the following parameters for the proposed method:σ=5.5, κ=1.4.Table 3
Region overlap metrics of different algorithms corresponding to Figure6 in the same order.
Image number
Algorithms
JS
DSC
FPR
FNR
Image #1
GAC
0.0033
0.0067
0.1095
0.9967
LSR
0.0026
0.0053
0.3212
0.9974
CV
0.0034
0.0069
0.1387
0.9966
FTC
0.0031
0.0061
0.1266
0.9845
LBF
0.0030
0.0059
0.1971
0.9970
LAC
0.0031
0.0062
0.1533
0.9969
LPFE
0.0023
0.0045
0.3869
0.9977
NFT
0.0012
0.0024
0.7226
0.9988
Our method
0.8919
0.9429
0.0365
0.0769
Image #2
GAC
0.0109
0.0215
0.0310
0.9891
LSR
0.0085
0.0169
0.0303
0.9915
CV
0.0090
0.0178
0.0207
0.9910
FTC
0.0119
0.0236
0.0205
0.9881
LBF
0.0099
0.0196
0.0455
0.9901
LAC
0.0117
0.0232
0.0038
0.9883
LPFE
0.0123
0.0242
0.0076
0.9877
NFT
0.0076
0.0151
0.0265
0.9924
Our method
0.8853
0.9392
0.0034
0.0573
Image #3
GAC
0.0973
0.1773
0.0807
0.9027
LSR
0.1055
0.1908
0.0088
0.8944
CV
0.1120
0.2014
0.3003
0.8880
FTC
0.0948
0.1732
0.4070
0.9052
LBF
0.1001
0.1820
0.0066
0.8998
LAC
0.0959
0.1750
0.0019
0.9041
LPFE
0.0914
0.1676
0.0016
0.9085
NFT
0.1124
0.2021
0.0082
0.8875
Our method
0.9580
0.9786
0.0010
0.0209In the next set of comparative experiments, we will give the convergence rate metrics for various methods, which are usually reflected by the number of iterations and the CPU time required to complete the entire evolution process. It is well known that the comparison of convergence rate metrics is meaningful only when the segmentation results of all the models participating in the comparison are correct. In view of this, we take the best approach when configuring constraints such as the initial level set function and the input parameters of the model. Here, we adopt the same comparison methods as the three sets of experiments mentioned above. Figure7 shows the results of our model and other comparison models using the same images and the same initial contours, where the first row shows the original images (they are a set of simple images with a uniform background and target) along with initial contours, and the second to the tenth rows are the segmentation results by using GAC, LSR, CV, FTC, LBF, LAC, LPFE, NFT and our method, respectively. From the final convergence curves, we can clearly see that all models output almost (there are only very small differences between the results) the same correct segmentation results. The corresponding convergence rate metrics for this set of experiments are shown in Table 4, and the sizes of these images are also shown under the image numbers. After analyzing the convergence rate metrics in Table 4, we find that our model achieves the best convergence performance because the graph cuts optimization method with very fast characteristics is adopted in its numerical solution process. Although FTC achieves only the second-best performance in the entire comparison methods family, it is the best among all the comparison models that use traditional numerical discretization methods; the reason is that it adopts a special element switching mechanism, which makes it unnecessary to solve the partial differential equation. The convergence performance is closely followed by other comparison models. Because they adopt more complex internal or external energy terms, they get a slower rate of convergence.Figure 7
The test images for convergence rate comparison experiments.Table 4
Comparison of convergence rate metrics for the images in Figure7 in the same order. The description format of the metrics is given as iteration number (CPU time (s)).
Image #1 (256 × 256)
Image #2 (256 × 256)
Image #3 (256 × 256)
GAC
300 (22.566)
270 (18.028)
259 (18.092)
LSR
59 (6.822)
150 (13.680)
53 (6.080)
CV
131 (7.432)
218 (14.258)
83 (5.518)
FTC
33 (3.826)
66 (3.181)
29 (3.399)
LBF
291 (19.130)
466 (31.517)
178 (12.007)
LAC
400 (51.588)
900 (97.469)
225 (28.699)
LPFE
159 (10.042)
191 (11.637)
139 (9.637)
NFT
394 (41.282)
407 (42.795)
316 (33.468)
Our method
9 (0.901)
18 (1.545)
8 (0.772)Through the above four sets of comparison experiments, we clearly see that the proposed model has the best performance in both segmentation accuracy and convergence rate under different image modalities, which is mainly attributed to the following factors:(1)
The proposed level set model possesses a strong ability of neighborhood modeling, thus it can effectively suppress the noise interference.(2)
This paper adopts an adaptive regularization strategy specifically for the zero level curve, thus this further enhances the noise suppression capability of the proposed model. At the same time, it also makes the proposed model has strong weak edge detection capability.(3)
Since graph cuts optimization method is used for numerical solution, the convergence rate of our model is greatly improved.


## 5. Conclusions

In this paper, we designed a robust MRF and adaptive regularization embedded level set segmentation model and solved it by graph cuts optimization. The existence of MRF makes the proposed method to have strong ability of neighborhood modeling, which makes it to have strong noise suppression performance. In addition, in order to further improve the noise suppression capability of the proposed model and accompany it with strong weak edge detection capability, we introduce a regularization strategy with adaptive characteristics. In numerical calculation stage, we use graph cuts optimization thought with extremely fast characteristics to realize the iterative solution of our model. In order to fully test the proposed model’s segmentation ability, we apply it to the segmentation tasks of different image modalities (including noisy synthetic images, medical images, infrared images, and clean synthetic images). Under the combined effect of the aforementioned three factors, the proposed model gets the best performance in both segmentation accuracy and convergence rate under different image modalities.



---

*Source: 8747385-2019-01-14.xml*