Magnetic Resonance Imaging Segmentation via Weighted Level Set Model Based on Local Kernel Metric and Spatial Constraint

Magnetic resonance imaging (MRI) segmentation is a fundamental and significant task since it can guide subsequent clinic diagnosis and treatment. However, images are often corrupted by defects such as low-contrast, noise, intensity inhomogeneity, and so on. Therefore, a weighted level set model (WLSM) is proposed in this study to segment inhomogeneous intensity MRI destroyed by noise and weak boundaries. First, in order to segment the intertwined regions of brain tissue accurately, a weighted neighborhood information measure scheme based on local multi information and kernel function is designed. Then, the membership function of fuzzy c-means clustering is used as the spatial constraint of level set model to overcome the sensitivity of level set to initialization, and the evolution of level set function can be adaptively changed according to different tissue information. Finally, the distance regularization term in level set function is replaced by a double potential function to ensure the stability of the energy function in the evolution process. Both real and synthetic MRI images can show the effectiveness and performance of WLSM. In addition, compared with several state-of-the-art models, segmentation accuracy and Jaccard similarity coefficient obtained by WLSM are increased by 0.0586, 0.0362 and 0.1087, 0.0703, respectively.


Introduction
Magnetic resonance imaging (MRI) is a part and parcel of medical imaging for its virtues such as rapid acquirement, non-intrusive and painless [1]. MRI images have been widely used in treatment evaluation, brain development monitoring, diagnosis, and so on [2][3][4][5]. However, the quality of MRI images is often influenced by various reasons such as low-contrast, noise and intensity inhomogeneity (IIH) during the imaging process, so it is not convenient to accurately segment and extract brain tissues. In addition, intensity inhomogeneity is called as bias field due to the property of slowly varying pixels in the same tissue [6]. Image segmentation has been extensively and deeply studied in computer vision due to its widespread application [7][8][9][10]. In the aspect of biomedical image analysis, it is a fundamental and complex task, which aims at assigning each pixel or voxel to the region with the same anatomical or biological meaning [11][12][13][14]. For example, by segmenting brain tissue, the doctors can detect and judge the changes of brain volume in the physiological or pathological state. Therefore, the need for methods that can accurately segment images has risen rapidly. Many segmentation methods have been presented by researchers to segment MRI images corrupted by intensity inhomogeneity, noise, and so on. These methods can be divided into several mainly categories: boundary-based methods [15][16][17][18][19], threshold-based methods [20][21][22], clustering methods [23][24][25][26][27], region growing methods [28,29], graph cuts methods [30][31][32] and deep learning method [33,34].
Deep learning is applied to medical image processing, segmentation and classification accuracy has been greatly improved, it effectively promotes the development of brain image automatic processing technology. However, this technology also has some problems.
Deep learning mostly depends on expert-labeled samples for training. Small-scale sample training has poor results. However, it is extremely difficult, expensive, and time-consuming to obtain large-scale expert-labeled brain segmentation samples. Most of the current studies have focused on the segmentation of glioma and meningioma, and there are few studies on segmentation for other brain lesions. Most deep learning networks are still in the experimental stage, and it will take a long time before clinical application. Because boundary-based methods do not need large-scale expert labeled samples, and the segmentation results have obvious and smooth contour. Therefore, boundary-based methods have been widely applied to MRI image segmentation.
Boundary-based methods depend on the intensity information of MRI image to detect the boundary of interested region [35]. Level set model is one of the most frequently used and well-established boundary-based methods for image segmentation and has attracted significant attention in the past few decades. It regards the contour called level set function as zero level set of higher dimensional function, which can transform the movement of contour into the evolution of level set function implicitly. Level set models can be subdivided into edge-based models and region-based models again [36].
Gradient information is employed by edge-based model as a constraint condition to detect the boundary of an interesting region. This method is particularly effective for the image with clear boundaries. However, for the image with weak boundaries, the level set function is usually difficult to converge to the correct boundaries due to the destruction of noise and intensity inhomogeneity [37]. Li et al. designed a distance regularized level set evolution (DRLSE) scheme [38], in which the image gradient information is used as driving force and regularization term is introduced into energy function to adjust the deviation between signed distance function and level set function to eliminate the defect of continuous re-initialization. Nevertheless, DRLSE only uses the image gradient information as a constraint condition to extract the boundary of each tissue. As a result, DRLSE is helpless to image with more than one target or weak boundaries. Then, Zhang et al. presented a level set evolution method driven by enhanced term to improve the effectiveness of DRLSE [39], the optimized area energy term was defined in their study to detect the boundaries of an image with several disjoint targets. Although the accuracy of this model is improved effectively, it is still sensitive to noisy pixels and initial contour.
Region-based models rely on the similarity of pixels belonging to the same anatomical tissue and relevant statistical information as the constraint conditions to detect the boundary of an interesting region. Region-based models have more advantages than edge-based models since they use region information rather than gradient information to drive the motion of level set function [40]. Therefore, they can be used to segment the images with weak boundaries and low contrast, and it is often insensitive to the position of the initial contour. The most representative model is Chan-Vese (CV) model [41]. It is not only suitable for the case of piecewise smooth but also for the case of piecewise constant. However, it presupposes that the intensity of each tissue is uniform, which means it cannot deal with the inhomogeneous MRI images well. Later, Li et al. a local clustering criterion in local intensity clustering (LIC) to process inhomogeneous intensity images and estimate bias field simultaneously [42]. LIC uses local intensity information and kernel function to homogenize the intensity of each tissue, which can effectively improve the accuracy of segmentation results. However, due to lack of spatial constraints, LIC is sensitive to noise, and the estimated bias field cannot satisfy the slowly varying property. Then, Feng et al. presented a local inhomogeneous intensity clustering (LINC) to improve the segmentation performance of LIC and ensure the properties of bias field [43]. More specifically, LINC is a LIC model with multiplicative intrinsic component optimization (MICO) [44], which utilizes clustering criterion defined in LIC and the basis function of bias field employed in MICO to estimate bias field. MICO employed the fourth-order Legendre polynomials as the basis function to model bias field. However, from the experimental results, LINC is still sensitive to noise, especially for images with weak boundaries. Recently, local information attracts more attention than global information since it can directly reflect the Entropy 2021, 23, 1196 3 of 20 similar property of the central pixel according to its neighbor pixel. Zhou et al. proposed a correntropy-based level set method (CLSM) for correcting bias field and segmenting medical images [45]. CLSM employed both local intensity information and bias field information to define a local bias-field-corrected image fitting energy, which is effective in segmenting inhomogeneous intensity medical images. Further, an optimal segmentation scheme was defined by Huang et al. using a fast level set model [46]. This model provided an optimal partition for every pixel in the entire image domain to inhomogeneous intensity images. In order to eliminate the influence of the bias field, MICO uses a fourth-order Legendre polynomial as the basis function to model bias field. In the process of energy minimization, it is a good modeling idea to find the best coefficient of the basis function to ensure the property of bias field. However, the usage of a fourth-order polynomial in its basis function will lead to higher computational complexity, and it is necessary to simplify the computational efficiency. Finally, the advantages and disadvantages of the relevant studied methods are summarized in Table 1. LIC [42] Able to estimate bias field and segment brain tissues simultaneously.
Bias field model is an idealized model without fully considering its own properties.
MICO [44] The slowly and the smoothly varying property of the bias field is ensured by a linear combination of a given set of smooth basis functions.
The model is not a level set method and is sensitive to noise without considering local neighborhood information.
LINC [43] A local clustering criterion function is defined to cluster intensities in the neighborhood for utilizing local neighborhood information.
All pixels including noise pixels are clustered into local clustering criterion, so LINC is sensitive to noise and weak boundaries.
CLSM [45] Incorporate the correntropy criterion into the energy function of local bias-field-corrected fitting image.
Difficult to discriminate pixels having same or minor differences between foreground and background, program execution efficiency is low.
Although many improvement measures have been proposed for the level set-based image segmentation scheme, there are still some shortcomings in the initialization of the level set function, the stability and accuracy of model evolution. Further, not only the intensity information but also more information such as variance, spatial position and graylevel should be incorporated into energy function to deal with MRI images. In this study, a weighted level set model (WLSM) based on local kernel metric and spatial constraint is proposed to segment brain MRI images corrupted by noise and intensity inhomogeneity. A neighborhood weighting method consisting of local variation, spatial distance and gray-level information is constructed to improve the resolution of low-contrast brain MRI images and reduce noise interference. Then, the coefficient of data term in energy function is enhanced by fusing fuzzy spatial constraint, which is more suitable for inhomogeneous MRI image segmentation. In addition, in order to ensure the stability of the level set function in the evolution process of level set, the single potential function used in the distance regularization term of the energy function is improved to a new double-well potential function. The organization of the rest of this study is as follows. Section 2 describes the WLSM proposed in this study in detail. Section 3 presents the experimental results consisting of visual experiments and quantitative evaluation. Finally, the conclusion of this study is given in Section 4.

Weighted Level Set Model
Brain MRI images are easily corrupted by noise and bias field, which makes it difficult for many existing algorithms to obtain satisfactory segmentation results. To this end, in response to practical problems, a weighted level set model combining local multipleinformation is proposed, which ensures accurate segmentation of brain MRI images by introducing kernel metrics and fuzzy spatial membership constraints.

Weighted Neighborhood Information
The improvement of image quality in this study mainly relies on the neighborhood information of the central pixel x i especially effective for the low-contrast brain MRI images, where white matter (WM) and gray matter (GM) are intertwined with each other and the boundary is blurred. Local variation, spatial distance and gray-level difference are used simultaneously to construct the weighted neighborhood information, which will be described separately in this subsection.

Local Variation Coefficient
The local variation depends on the intensity mean and variance of pixels fall into the local window of central pixel x i . Variance, which is defined as the degree of deviation between pixels and mean, can reflect the difference between pixels intuitively. Therefore, the local variation coefficient C ij of each neighbor pixel x j for the center pixel x i is defined as: where x j denotes the gray scale of neighborhood pixel of x j , that is, x j is the neighbor pixel of x i , while x j is the neighbor pixel of x j . var(x) and x stand for the variance and mean of the corresponding image block, respectively. N j denotes the neighborhood centered on x j . The greater the difference between pixels x j , the larger the value of C ij , which means that it reflects the degree of homogeneity of pixels x j . If x j is similar to x j , such as the region with homogeneous intensity, the value of C ij will be small, and vice versa. Thus, it can utilize more contextual information to explore the local variation because it is computed in the neighborhood of x j . Then, C ij is projected into the kernel metric space, and its larger value will result in a smaller gray level weight because of the fast damping property of kernel function. Both Gaussian kernel function and Laplace kernel function have the ability to restrain noise, but the performance of Gaussian kernel function depends largely on the choice of coefficient. Laplace kernel function is the variant of Gaussian kernel function, and it is insensitive to the choice of coefficient. Therefore, the weighting coefficient γ ij of each pixel x j in the neighborhood of x i using Laplace kernel can be expressed as: where l is a bandwidth constant. N i is the neighborhood centered on x i .

Adaptive Spatial Measure
The local spatial measure used in some typical studies is still the Euclidean distance due to its computationally simple. However, Euclidean distance will easily result in incorrect measure results for the images destroyed by noise, intensity inhomogeneity, and so on. Therefore, a nonlinear version of the linear distance is constructed using kernel metric to improve Euclidean distance, which has been used in fuzzy spatial C-means [47] and support vector clustering [48]. Gaussian kernel function is an efficient distance measure, which is given by: where σ is the bandwidth coefficient that is often set to a constant in most studies. According to the above description about Gaussian kernel function, it is sensitive to the choice of coefficient and σ has a significant effect on its performance, so σ should not be fixed as a constant. Thus, WLSM enhances the Gaussian kernel function by calculating σ based on the average variance of the distance between x j and x i to adaptively update the value of σ. Let Euclidean distance between x j and x i is: The degree of deviation between x j and x i can be seen from d ij . Then, the mean distance d of d ij is given by: where n is the number of x j . Accordingly, σ i can be calculated based on mean square variance of these distances Finally, the enhanced Gaussian kernel distance is rewritten as: Comparing Equations (3) and (7), the constant σ is replaced by the variable σ i . In order to visually display the role of D ij , take an eight-neighborhood image patch in the white matter region of brain MRI image, and the value of D ij after normalization is shown in Figure 1. There are two pixels in the image patch are destroyed in Figure 1a, and the weights of noisy pixel and the pixel with the same intensity obtained from enhanced Gaussian kernel distance in Figure 1b are greatly different. This is partly because the damping degree of kernel distance will decrease the weight of noisy pixels, but mainly because the bandwidth coefficient σ i of the enhanced Gaussian kernel distance can be calculated adaptively on the basis of the mean square variance of the distance between neighbor pixels and the central pixel. the damping degree of kernel distance will decrease the weight of noisy pixels, but mainly because the bandwidth coefficient σi of the enhanced Gaussian kernel distance can be calculated adaptively on the basis of the mean square variance of the distance between neighbor pixels and the central pixel.

Synthetic Weight of Neighborhood Term
The local gray level information x j with gray level weight γ ij and Gaussian kernel function D 2 ij can be used to update the current central pixel x i .
An original image with 5% noise obtained from BrainWeb [49] and its corresponding enhanced image using weighted neighborhood information are compared in Figure 2, and four partial enlarged areas are also displayed. Due to the noise in the original image, the edges of WM and GM are intertwined, it is difficult to distinguish the edge of different brain tissue. However, the enhanced image is much more noise free than the original image. Further, the partially enlarged region can show more details of the boundaries of WM and GM. Therefore, WLSM can improve the quality of images effectively and reduce the influence of noisy pixels.

Synthetic Weight of Neighborhood Term
The local gray level information xj with gray level weight γij and Gaussian kernel function can be used to update the current central pixel xi.
An original image with 5% noise obtained from BrainWeb [49] and its corresponding enhanced image using weighted neighborhood information are compared in Figure  2, and four partial enlarged areas are also displayed. Due to the noise in the original image, the edges of WM and GM are intertwined, it is difficult to distinguish the edge of different brain tissue. However, the enhanced image is much more noise free than the original image. Further, the partially enlarged region can show more details of the boundaries of WM and GM. Therefore, WLSM can improve the quality of images effectively and reduce the influence of noisy pixels.

The Improved External Energy Function
An image I(x) is divided into N disjoint regions Ω1, …, ΩN, and the local intensity clustering defined in our previous study [35] described that pixels x around y with radius ρ is given by Oy = {x: |x -y| ≤ ρ}. Thus, all pixels in the neighborhood Oy can be divided into k classes, which is given by: In other words, this classification criterion means that intensities fall into Oy can be classified into N clusters. Subsequently, the clustering function defined based on such

The Improved External Energy Function
An image I(x) is divided into N disjoint regions Ω 1 , . . . , Ω N , and the local intensity clustering defined in our previous study [35] described that pixels x around y with radius ρ is given by O y = {x: |x − y| ≤ ρ}. Thus, all pixels in the neighborhood O y can be divided into k classes, which is given by: In other words, this classification criterion means that intensities fall into O y can be classified into N clusters. Subsequently, the clustering function defined based on such criterion is also suitable for this study to classify intensities. According to Ref. [38], the date term ε constituting the external energy function of the level set can be expressed as: where K(y − x) is a function used to control the range of O y , which has been defined in LIC and should satisfy K(s)ds = 1. c k is the k-th clustering center. w T G(y) denotes the bias field b(x), where G(y) is the combination of Legendre polynomial functions that can ensure b(x) is slowly changing and vector w = (ω 1 , . . . , ω N ) T is the optimal coefficient. I(x) denotes the original image, but the weighted neighborhood information can improve the quality of image before the energy function is minimized, so the clustering function can be re-represented as: where I(x) denotes the enhanced image of I(x) using Equation (8). In Ref. [39], a level set function φ(x) can represent two disjoint regions Ω 1 and Ω 2 by A multi-phase level set function Φ = (φ 1 , . . . , φ i ) can denote N disjoint regions Ω 1 , . . . , Ω N by corresponding membership functions M k (Φ). For example, φ 1 (x) and ). Therefore, clustering Function (11) can be converted to the data term of the level set model as follows: Clustering center c k is replaced by the vector c = (c 1 , . . . , c N ) for convenience. Therefore, the data term ε can be re-represented as a new expression ε(Φ, w, c) about level set function Φ, optimal coefficient w and clustering center vector c. Swapping the order of integrations will not change the data term, so it is written as: where λ k is the weighted coefficient and e k is expressed as: Its numerical expression is given by: where 1 K denotes K(y − x)dy. The weighted coefficient λ k of data term is often set as λ 1 = . . . = λ k = 1 in some typical level set models, such as LIC and LINC. Nevertheless, without any constraints, Φ cannot reach or directly pass through the weak boundaries corrupted by low-contrast and thus leads to the interesting region which cannot be completely extracted. Therefore, λ k should not be set as a global constant at least in the vicinity of zero level set. According to AFLSM [30], λ k has two significant roles, one is the direction of Φ depends on the value of λ k , the other is to control the speed of Φ. In other words, Φ should be accelerated when it is far away from the boundaries of the interested region and vice versa. Therefore, λ k should be improved via contextual information. Fuzzy spatial constraint obtained by fuzzy clustering is suitable for improving λ k because it assumes pixels in the vicinity of boundaries belong to different clusters according to their membership matrix u ki , where k denotes the k-th clustering center and i is the i-th pixel. Then, the modified weighted coefficient λ k constrained by u ki is defined as: The difference between the λ k defined in this study and the one defined in our previous study AFLSM is mainly reflected in two aspects. For one thing, the coefficient α of λ k defined in AFLSM needs to be fine-tuned based on images with different characteristics, but the λ k defined in WLSM completely depends on the membership matrix u ki , which means the modified λ k in this study is fully adaptive and does not need to be adjusted manually. For another, all λ k defined in WLSM are normalized to ensure the stability of numerical calculation and control the variation of level set function exactly. Thus, the modified λ k is superior to the one defined in AFLSM and the segmentation accuracy of WLSM can be further improved for the images with weak boundaries.

Internal Energy Function
Li et al. regularized the level set function in LBF [50] to a smooth zero level set by computing the length of its contour, which is defined as the following length term: where ∇ is the gradient operator. A conventional level set function has to be reinitialized during each iteration, which is a time-consuming process, especially for a model with thousands of iterations. Accordingly, LBF presented a distance regularization term P(Φ) that utilizes the unique property |∇Φ| = 1 of signed distance functions to avoid the expensive re-initialization during the process of evolution. That is, the distance regularization term can drive Φ to approach the signed distance function to eliminate the deviation, which is defined by: where p(s) is a potential function. p(s) is often set as p(s) = (s − 1) 2 /2, where s = 1 is its minimum point used to maintain the property of signed distance functions. However, the single potential function is not always effective in practical applications because of its forwardand-backward diffusion and thus cannot steadily maintain the property |∇Φ| = 1 [33]. Therefore, WLSM improved the single potential function by introducing two minimum points s = 0 and s = 1, which is aimed to solve the drawbacks of LBF and maintain the property |∇Φ| = 1 only in a vicinity of zero level set contour. Such a potential function is called the double-well potential function p 2 (s) because of its two minimum points, which should satisfy the following four conditions: (1) Two minimum points of double potential function p 2 (s) should be at s = 0 and s = 1, respectively; Based on the above four constraints, p 2 (s) is defined in WLSM can be expressed as: The curve of p 2 (s) and d p (s) are shown in Figure 3. When p 2 (s) = 0, it can be seen from Figure 3a that p 2 (s) have two minimum points (no other minimum points) at s = 0 and s = 1, respectively. The first derivative p 2 (s) of p 2 (s) is given by: The second derivative p 2 "(s) of p 2 (s) is given by: Entropy 2021, 23, 1196 9 of 20 which means the novel p 2 (s) is second-order derivable in [0,∞). Next, the function d p (s) based on p 2 (s) is defined as: The curve of p2(s) and dp(s) are shown in Figure 3. When p2(s) = 0, it can be seen from Figure 3a that p2(s) have two minimum points (no other minimum points) at s = 0 and s = 1, respectively. The first derivative p2′(s) of p2(s) is given by: The second derivative p2″(s) of p2(s) is given by: which means the novel p2(s) is second-order derivable in [0,∞). Next, the function dp(s) based on p2′(s) is defined as: As shown in Figure 3b, it is not difficult to prove that dp(s) derived from p2(s) also satisfies the above conditions Equations (3) and (4), which are to ensure the diffusion rate of p2(s) is bounded. Further, the novel double-well potential function p2(s) not only can eliminate the drawbacks of single potential function but also can penalize its deviation from the signed distance functions. Therefore, WLSM with double potential function p2(s) is more suitable for segmenting images with noise and intensity inhomogeneity.

Energy Formulation and Its Minimization
According to the external function and the internal function, the entire energy formulation can be expressed as: where ε (Φ, w, c) is the data term defined in Equation (14), ν and μ are the weighted coefficients of length term L(Φ) defined in Equation (18) and distance regularization term P(Φ) defined in Equation (19), respectively. The method to minimize Φ is to solve the following gradient descent flow equation As shown in Figure 3b, it is not difficult to prove that d p (s) derived from p 2 (s) also satisfies the above conditions Equations (3) and (4), which are to ensure the diffusion rate of p 2 (s) is bounded. Further, the novel double-well potential function p 2 (s) not only can eliminate the drawbacks of single potential function but also can penalize its deviation from the signed distance functions. Therefore, WLSM with double potential function p 2 (s) is more suitable for segmenting images with noise and intensity inhomogeneity.

Energy Formulation and Its Minimization
According to the external function and the internal function, the entire energy formulation can be expressed as: where ε(Φ, w, c) is the data term defined in Equation (14), ν and µ are the weighted coefficients of length term L(Φ) defined in Equation (18) and distance regularization term P(Φ) defined in Equation (19), respectively. The method to minimize Φ is to solve the following gradient descent flow equation where ∂F/∂Φ is Gâteaux derivative and t is a time variable. The right-hand side of the equation shows the steepest descent direction of the energy function F. Accordingly, the minimization of Φ is given by: where e k has been given in Equation (16). The Gâteaux derivative of length term L(Φ) is given by: where δ(φ) is the derivative of H(φ) and div(·) is divergence operator. Then, the Gâteaux derivative of distance regularization term P(Φ) is given by: where the function d p (s) has been defined in Equation (23). The method to minimize optimal coefficient w is to solve the following equation: where v is a column vector defined by: and A is a matrix defined by: Thus, according to the expression of v and A, the minimization of w is given by: The minimization of c k is given by: The implementation process of WLSM can be expressed as Algorithm 1.

Experiments and Results
In this section, WLSM is used to segment synthetic and real MRI images and compared with state-of-the-art models. Firstly, it is used to segment MRI images with varying degrees of noise to demonstrate that the weighted neighborhood information embedded in WLSM can enhance the quality of MRI images and improve segmentation accuracy (SA). Then, it is used to correct inhomogeneous intensity MRI images and estimate bias field. Finally, the sagittal, coronal and axial slices of synthetic and real MRI images are segmented by WLSM to show its performance from different aspects. BrainWeb is a simulated brain database (SBD) and includes a set of data generated by an MRI simulator that is close to real brain MR images [49], and it was developed by the brain imaging center of the Montreal Neurology Institute at McGill University. The brain MR image data in SBD is composed of a three-dimensional matrix of 181 × 217 × 181 voxels, which can simulate T1, T2 and PD (proton-density) weighted brain MRI images, and the slice thickness, noise level and intensity inhomogeneity can be set by oneself. According to the requirement of the research task, 3D images can be sliced and extracted from three planes (sagittal, coronal and axial) to obtain 2D image data. The real images used in this study are from the internet brain segmentation repository (IBSR) database provided by Formal Measurement Center (CMA) of Massachusetts General Hospital [51]. IBSR database is a set of clinical data generated by real MRI scans, which contains different levels of noise, and there are also varying degrees of intensity inhomogeneity, covering various problems that may arise in real MR data segmentation. Unless otherwise specified, parameters used in WLSM are fixed as ν = 0.001 × 255 × 255, µ = 1, and ∆t = 0.1. All experimental results were implemented by Matlab R2019a on a computer with Intel Core i5-8300H 2.3 GHz CPU, 8 GB RAM, Windows 10 operating system, and the segmentation results with final zero-level contours of φ 1 in red and φ 2 in blue are shown.

Results on Noisy Images
First of all, WLSM model is applied to 20 T1-weighted normal brain MRI images with 181 × 217 pixels, 1 mm slice thickness, 5% noise, without intensity inhomogeneity, their cranium and blood vessels have been dislodged before segmentation processing. In Figure 4, serial numbers of the sliced images from the axial plane are 75 and 85, respectively. Original images in the first column are destroyed by noise especially in the regions that WM and GM are interweaved with each other, which can be easily seen from the partial enlarged regions marked by the green rectangles. For comparison, the enhanced images using weighted neighborhood information in the second column are more noise free than original images and the boundary of each tissue can be clearly distinguished because of the constraints of weighted neighborhood information. To be specific, WLSM utilizes the neighbor information consisting of local variation, spatial distance and gray-level information to update the current central pixel. The neighbor pixels in the same tissue will occupy bigger weights and vice versa. Therefore, pixels located in the weak boundaries can be categorized into the corresponding tissue according to their weighted neighborhood information. Subsequently, benefitting from such the weighted neighborhood information, level set functions are able to accurately extract the boundary of brain tissue, which can be observed from the third column, the last column is the ground truth.  The higher the noise level, the more difficult it is to accurately segment the noisy MRI images. Thus, brain MRI images with the difference of noise level are segmented by WLSM to testify its robustness to noise. The sliced images with 3%, 5%, 7% and 9% noise are displayed in the first to the fourth columns of Figure 5, respectively. The quality of original images in the first row gradually decreases with the increase of noise, which will easily result in non-robust segmentation results. However, the segmentation results obtained by WLSM in the second row have no obvious discrepancy even in the case of severe noise, which can show its sturdy robustness to the different noise levels. The higher the noise level, the more difficult it is to accurately segment the noisy MRI images. Thus, brain MRI images with the difference of noise level are segmented by WLSM to testify its robustness to noise. The sliced images with 3%, 5%, 7% and 9% noise are displayed in the first to the fourth columns of Figure 5, respectively. The quality of original images in the first row gradually decreases with the increase of noise, which will easily result in non-robust segmentation results. However, the segmentation results obtained by WLSM in the second row have no obvious discrepancy even in the case of severe noise, which can show its sturdy robustness to the different noise levels.
(a) (b) (c) (d) The higher the noise level, the more difficult it is to accurately segment the noisy MRI images. Thus, brain MRI images with the difference of noise level are segmented by WLSM to testify its robustness to noise. The sliced images with 3%, 5%, 7% and 9% noise are displayed in the first to the fourth columns of Figure 5, respectively. The quality of original images in the first row gradually decreases with the increase of noise, which will easily result in non-robust segmentation results. However, the segmentation results obtained by WLSM in the second row have no obvious discrepancy even in the case of severe noise, which can show its sturdy robustness to the different noise levels. Because fuzzy C-means (FCM) clustering and its variants are often used to segment brain MRI images, and in the WLSM model, the spatial neighborhood constraint is also realized by the idea of the FCM. Therefore, to further compare the segmentation performance, the two clustering algorithms are also compared. One is the standard FCM algorithm, and the other is the adaptively regularized kernel-based fuzzy C-means clustering (ARKFCM) algorithm in Ref. [52]. Figure 6 shows the comparison results of FCM, LINC [38], MICO [39], ARKFCM [52] LIC [37] and WLSM on brain MRI images with 5% noise. FCM is very sensitive to noise, and the segmentation effect is poor. LINC employs the local clustering function to define the energy function, which can cluster inhomoge- Because fuzzy C-means (FCM) clustering and its variants are often used to segment brain MRI images, and in the WLSM model, the spatial neighborhood constraint is also realized by the idea of the FCM. Therefore, to further compare the segmentation performance, the two clustering algorithms are also compared. One is the standard FCM algorithm, and the other is the adaptively regularized kernel-based fuzzy C-means clustering (ARK-FCM) algorithm in Ref. [52]. Figure 6 shows the comparison results of FCM, LINC [38], MICO [39], ARKFCM [52] LIC [37] and WLSM on brain MRI images with 5% noise. FCM is very sensitive to noise, and the segmentation effect is poor. LINC employs the local clustering function to define the energy function, which can cluster inhomogeneous intensity in the neighborhood. However, all pixels including noisy pixels are clustered in the local clustering function, so the performance of LINC is easily influenced by noise. Accordingly, it can be seen from Figure 6c that LINC obtained incorrect results in the face of images corrupted by noise. MICO is based on the hard clustering to segment images, where hard clustering assumes that each pixel in the entire image domain only belongs to one cluster, so it is also sensitive to noise. Thus, the segmentation results obtained by MICO in Figure 6d are easily influenced by noise pixels. In addition, the results have numerous wrong contours representing noisy pixels. ARKFCM employs the heterogeneity of grayscales in the neighborhood and exploit this measure for local contextual information and replace the standard Euclidean distance with Gaussian radial basis kernel functions, which can effectively suppress the noise in the image, but ARKFCM cannot estimate the bias field in brain MRI images, and it is susceptible to the interference of intensity inhomogeneity, as shown in Figure 6e. LIC uses the intensity clustering property to define a global criterion function, which can estimate bias field and segment images simultaneously. However, it is also sensitive to noise because of the lack of spatial constraints. As shown in Figure 6f, LIC cannot accurately extract the boundary between GM and WM and thus leads to over-segmentation. By contrast, the segmentation results by WLSM displayed in Figure 6g shows superior performance than the other five models.
found by the algorithm and Gi is the pixel number belonging to the ith cluster in GT. The higher the value of SA means that model can obtain more accurate results on the segmented images. In this experiment, we selected 10 groups of brain MRI images as experimental samples to compare the segmentation results. SA of FCM, LINC, MICO, ARKFCM, LIC and WLSM are shown in Table 2. Both the SA of WM and GM and the average SA can show the superiority of WLSM, which objectively proves that WLSM obtains more accurate segmentation results on noisy MRI images.   In order to quantitatively analyze the experimental results, segmentation accuracy (SA) is often used to evaluate the segmentation performance of algorithms, which is defined as:

Results on Inhomogeneous Intensity Images
where k is the number of clusters, S i is the pixel number belonging to the ith cluster found by the algorithm and G i is the pixel number belonging to the ith cluster in GT. The higher the value of SA means that model can obtain more accurate results on the segmented images. In this experiment, we selected 10 groups of brain MRI images as experimental samples to compare the segmentation results. SA of FCM, LINC, MICO, ARKFCM, LIC and WLSM are shown in Table 2. Both the SA of WM and GM and the average SA can show the superiority of WLSM, which objectively proves that WLSM obtains more accurate segmentation results on noisy MRI images.

Results on Inhomogeneous Intensity Images
In general, MRI images are often corrupted by noise and intensity inhomogeneity simultaneously during the process of imaging, so the ability to correct inhomogeneous intensity and estimating bias field is the part and parcel of the performance of the model. Therefore, in the experiment of this subsection, WLSM is used to segment brain MRI images with 3% noise and varying degrees of intensity inhomogeneity to testify its performance. Images with 3% noise and 60%, 80% and 100% intensity inhomogeneity are displayed on the first row to the third row of Figure 7, respectively. The corrected images in the second column are much more homogeneous than one of the original images in the first column, and the segmentation results in the third column have no obvious discrepancy whether in the case of lower inhomogeneous intensity or higher inhomogeneous intensity, which illustrates that WLSM can correct inhomogeneous intensity and segment images simultaneously without being influenced by noise and intensity inhomogeneity. Benefitting from the third-order orthogonal Legendre functions, the bias field in the fourth column estimated by WLSM satisfies the property of smooth changing. Then, the histograms of original and corrected images are plotted in Figure 8 to objectively compare the images quality. The histograms of original images have no obvious peaks that represent corresponding tissues because of the influence of noise and intensity inhomogeneity. In contrast, there are well-defined peaks in histograms of the corrected images, which objectively reflect the improvement of the image quality. mogeneity. Benefitting from the third-order orthogonal Legendre functions, the bias field in the fourth column estimated by WLSM satisfies the property of smooth changing. Then, the histograms of original and corrected images are plotted in Figure 8 to objectively compare the images quality. The histograms of original images have no obvious peaks that represent corresponding tissues because of the influence of noise and intensity inhomogeneity. In contrast, there are well-defined peaks in histograms of the corrected images, which objectively reflect the improvement of the image quality.  In this experiment, because FCM and ARKFCM algorithms do not have the ability to estimate the bias field, we only compare the segmentation results of LINC, MICO, LIC and WLSM. In Figure 9, the tested brain MRI images are corrupted by 5% noise and 100% intensity inhomogeneity. Original images, corrected images, segmentation results and bias field obtained by each model are shown in the first to the fourth rows, respectively. As shown in the first and second columns of the second row, the corrected images obtained by LINC and MICO are still inhomogeneous especially in the regions that WM and GM are interweaved with each other. It can be seen from the third column of the last row that LIC cannot ensure the property of the bias field due to a lack of smooth basis functions. The corrected image obtained by WLSM is more homogeneous than the original image and the one obtained by LINC, MICO and LIC, respectively. In addition, the third-order orthogonal Legendre functions used by WLSM can ensure the especial property of bias field. Therefore, WLSM can estimate bias field excellently and the cor- In this experiment, because FCM and ARKFCM algorithms do not have the ability to estimate the bias field, we only compare the segmentation results of LINC, MICO, LIC and WLSM. In Figure 9, the tested brain MRI images are corrupted by 5% noise and 100% intensity inhomogeneity. Original images, corrected images, segmentation results and bias field obtained by each model are shown in the first to the fourth rows, respectively. As shown in the first and second columns of the second row, the corrected images obtained by LINC and MICO are still inhomogeneous especially in the regions that WM and GM are interweaved with each other. It can be seen from the third column of the last row that LIC cannot ensure the property of the bias field due to a lack of smooth basis functions.
The corrected image obtained by WLSM is more homogeneous than the original image and the one obtained by LINC, MICO and LIC, respectively. In addition, the third-order orthogonal Legendre functions used by WLSM can ensure the especial property of bias field. Therefore, WLSM can estimate bias field excellently and the corrected image will much more homogeneous after removing the estimated bias field. Then, the iterations and calculative time of LINC, MICO, LIC and WLSM are contrasted and the results are shown in Table 3. The size of the original image in the first column of Figure 9 is 181 × 217. It can be observed from the last row of Table 3 that WLSM consumes less calculative time than the other three models in segmenting brain MRI images corrupted by both noise and intensity inhomogeneity.  Jaccard similarity coefficient (JSC), the measure of similarity between two sets, is applied to quantitatively compare the performance of segmentation model. The definition of JSC is given by: where S is the segmentation result obtained by each model and SG is GT. The higher value of JSC indicates that the corresponding model is better than other methods. To quantitatively compare the accuracy of LINC, MICO, LIC and WLSM, the JSC obtained by the above four models is plotted in the box plots of Figure 10. As can be seen, the JSC of WM and GM obtained by WLSM exhibit the higher values among the four models, which can quantitatively demonstrate that WLSM can reach higher segmentation accuracy in comparison with LINC, MICO and LIC.  Jaccard similarity coefficient (JSC), the measure of similarity between two sets, is applied to quantitatively compare the performance of segmentation model. The definition of JSC is given by: where S is the segmentation result obtained by each model and S G is GT. The higher value of JSC indicates that the corresponding model is better than other methods. To quantitatively compare the accuracy of LINC, MICO, LIC and WLSM, the JSC obtained by the above four models is plotted in the box plots of Figure 10. As can be seen, the JSC of WM and GM obtained by WLSM exhibit the higher values among the four models, which can quantitatively demonstrate that WLSM can reach higher segmentation accuracy in comparison with LINC, MICO and LIC.

Results on Sagital, Coronal and Axial Slices of Images
In this subsection, the sagittal, coronal and axial slices of synthetic and real brain MRI images are segmented to show the effectiveness of models. Firstly, the segmentation results of synthetic brain MRI images obtained by FCM, LINC, MICO, ARKFCM, LIC and WLSM are shown in Figure 11. The sliced images from the sagittal, coronal and axial planes are displayed in the first to the third rows of Figure 11, respectively. Original images in the first column are corrupted by noise and intensity inhomogeneity simultaneously. FCM is not robust to noise and bias field, LINC still cannot accurately extract the interested tissues and thus leads to segmentation errors. MICO cannot distinguish the boundary of WM and GM, because their gray values are too similar in the case of severe inhomogeneous intensity, which can be seen from the sagittal and axial slice images in the second and third rows of the fourth column. ARKFCM cannot eliminate the influence of the bias field in the brain MRI images, and it has encountered a problem in the segmentation of white matter, as shown in the fifth column. LIC is sensitive to noise pixels because of the deficiency of spatial constraints. Hence, the segmentation results shown in the sixth column are still corrupted by noise especially in the boundary of WM and GM. As shown in the seventh column, WLSM can accurately extract the boundary of WM and GM whether on the sagittal, coronal or corrupted slice images, which means WLSM is effective on synthetic images with different planes.

Results on Sagital, Coronal and Axial Slices of Images
In this subsection, the sagittal, coronal and axial slices of synthetic and real brain MRI images are segmented to show the effectiveness of models. Firstly, the segmentation results of synthetic brain MRI images obtained by FCM, LINC, MICO, ARKFCM, LIC and WLSM are shown in Figure 11. The sliced images from the sagittal, coronal and axial planes are displayed in the first to the third rows of Figure 11, respectively. Original images in the first column are corrupted by noise and intensity inhomogeneity simultaneously. FCM is not robust to noise and bias field, LINC still cannot accurately extract the interested tissues and thus leads to segmentation errors. MICO cannot distinguish the boundary of WM and GM, because their gray values are too similar in the case of severe inhomogeneous intensity, which can be seen from the sagittal and axial slice images in the second and third rows of the fourth column. ARKFCM cannot eliminate the influence of the bias field in the brain MRI images, and it has encountered a problem in the segmentation of white matter, as shown in the fifth column. LIC is sensitive to noise pixels because of the deficiency of spatial constraints. Hence, the segmentation results shown in the sixth column are still corrupted by noise especially in the boundary of WM and GM. As shown in the seventh column, WLSM can accurately extract the boundary of WM and GM whether on the sagittal, coronal or corrupted slice images, which means WLSM is effective on synthetic images with different planes.

Results on Sagital, Coronal and Axial Slices of Images
In this subsection, the sagittal, coronal and axial slices of synthetic and real brain MRI images are segmented to show the effectiveness of models. Firstly, the segmentation results of synthetic brain MRI images obtained by FCM, LINC, MICO, ARKFCM, LIC and WLSM are shown in Figure 11. The sliced images from the sagittal, coronal and axial planes are displayed in the first to the third rows of Figure 11, respectively. Original images in the first column are corrupted by noise and intensity inhomogeneity simultaneously. FCM is not robust to noise and bias field, LINC still cannot accurately extract the interested tissues and thus leads to segmentation errors. MICO cannot distinguish the boundary of WM and GM, because their gray values are too similar in the case of severe inhomogeneous intensity, which can be seen from the sagittal and axial slice images in the second and third rows of the fourth column. ARKFCM cannot eliminate the influence of the bias field in the brain MRI images, and it has encountered a problem in the segmentation of white matter, as shown in the fifth column. LIC is sensitive to noise pixels because of the deficiency of spatial constraints. Hence, the segmentation results shown in the sixth column are still corrupted by noise especially in the boundary of WM and GM. As shown in the seventh column, WLSM can accurately extract the boundary of WM and GM whether on the sagittal, coronal or corrupted slice images, which means WLSM is effective on synthetic images with different planes.  Finally, WLSM is applied to segment the sagittal, coronal and axial slices of real brain MRI images and is also compared with FCM, LINC, MICO, ARKFCM and LIC to further demonstrate the effectiveness of WLSM. The segmentation results of real images obtained from the four models are shown in Figure 12. The real sagittal, coronal and axial slices are shown in the first to the third rows of Figure 12, respectively. Compared with all the experimental results, WLSM is more accurate than FCM, LINC, MICO, ARKFCM and LIC, which illustrates that it can be employed to segment real brain MRI images with different planes without being influenced by the skull. Finally, WLSM is applied to segment the sagittal, coronal and axial slices of real brain MRI images and is also compared with FCM, LINC, MICO, ARKFCM and LIC to further demonstrate the effectiveness of WLSM. The segmentation results of real images obtained from the four models are shown in Figure 12. The real sagittal, coronal and axial slices are shown in the first to the third rows of Figure 12, respectively. Compared with all the experimental results, WLSM is more accurate than FCM, LINC, MICO, ARKFCM and LIC, which illustrates that it can be employed to segment real brain MRI images with different planes without being influenced by the skull.

Discussion
The edge-based level set model is to segment and extract the corresponding brain tissue according to the boundary of the target. Therefore, for images with weak tissue boundaries or low contrast, the conventional level set function cannot accurately identify the boundaries of each tissue. Firstly, the WLSM algorithm targets the weak boundary of low contrast images and proposes a method to improve the edge contrast of brain tissue. This method uses local variance, spatial distance and gray information to update the current central pixel. The updated image has high contrast at the weak boundary of the brain tissues, so it lays a sound foundation for the subsequent level set function to extract the target edge. Secondly, the third-order Legendre polynomial function with orthogonality is used as the basis function to estimate the bias field, and its optimal coefficient is found in the iterative process. Compared with the traditional linear bias field estimation model, WLSM can ensure that bias field variation is smooth and slow. Thirdly, taking the membership function of FCM as the adaptive coefficient of the data item in the energy function changes the defect that the original coefficient is constant. It can not only overcome the problem that level set function is sensitive to initialization but also adaptively control the evolution of level set function. Finally, the single potential function in distance regularization term is replaced by a novel double-well potential function, which can effectively avoid reinitialization and maintain the accuracy and stability of the evolution of the level set model.
Of course, the proposed algorithm still needs to be improved. In the application of image segmentation, WLSM model only segments normal brain MRI images, and does not process brain MRI images with diseases such as multiple sclerosis, brain tumors and further research can be carried out on such images. Furthermore, this paper only studies the segmentation algorithm of 2D MRI slice images and does not directly segment 3D MRI images. In the future, we will further improve WLSN model combining deep

Discussion
The edge-based level set model is to segment and extract the corresponding brain tissue according to the boundary of the target. Therefore, for images with weak tissue boundaries or low contrast, the conventional level set function cannot accurately identify the boundaries of each tissue. Firstly, the WLSM algorithm targets the weak boundary of low contrast images and proposes a method to improve the edge contrast of brain tissue. This method uses local variance, spatial distance and gray information to update the current central pixel. The updated image has high contrast at the weak boundary of the brain tissues, so it lays a sound foundation for the subsequent level set function to extract the target edge. Secondly, the third-order Legendre polynomial function with orthogonality is used as the basis function to estimate the bias field, and its optimal coefficient is found in the iterative process. Compared with the traditional linear bias field estimation model, WLSM can ensure that bias field variation is smooth and slow. Thirdly, taking the membership function of FCM as the adaptive coefficient of the data item in the energy function changes the defect that the original coefficient is constant. It can not only overcome the problem that level set function is sensitive to initialization but also adaptively control the evolution of level set function. Finally, the single potential function in distance regularization term is replaced by a novel double-well potential function, which can effectively avoid reinitialization and maintain the accuracy and stability of the evolution of the level set model.
Of course, the proposed algorithm still needs to be improved. In the application of image segmentation, WLSM model only segments normal brain MRI images, and does not process brain MRI images with diseases such as multiple sclerosis, brain tumors and further research can be carried out on such images. Furthermore, this paper only studies the segmentation algorithm of 2D MRI slice images and does not directly segment 3D MRI images. In the future, we will further improve WLSN model combining deep learning to segment images in 3D and try to apply the multiphase formulation to segment brain tissues on public image repositories.

Conclusions
In this study, according to the deficiencies in the initialization of level set function, the stability and accuracy of model evolution in the current segmentation model, a weighted level set model (WLSM) based on local kernel metric and spatial constraints is proposed to segment brain MRI images corrupted by noise and intensity inhomogeneity. WLSM can not only accurately control the evolution of level set function, but also overcome the problem of model reinitialization. Compared with state-of-the-art models, the visual experiments including both synthetic and real MRI images demonstrate the superiority of WLSM.