WANG Yan,XIANG Yongjia,ZHANG Xuyuan,and WU Dan

1 College of Biomedical Engineering & Instrument Science, Zhejiang University,Hangzhou 310027,China.

2 School of Mathematical Sciences,Chongqing Normal University,Chongqing 401331,China.

3 College of Electronics and Information Engineering,Shenzhen University,Shenzhen 518060,China.

Abstract. Segmentation of images with intensity inhomogeneity is a significant task in the field of image processing,especially in medical image processing and analysis.Some local region-based models work well on handling intensity inhomogeneity, but they are always sensitive to contour initialization and high noise. In this paper, we present an adaptive segmentation model for images with intensity inhomogeneity in the form of partial differential equation. Firstly, a global intensity fitting term and a local intensity fitting term are constructed by employing the global and local image information,respectively. Secondly,a tradeoff function is defined to adjust adaptively the weight between two fitting terms,which is based on the neighborhood contrast of image pixel. Finally,a weighted regularization term related to local entropy is used to ensure the smoothness of evolution curve. Meanwhile,a distance regularization term is added for stable level set evolution. Experimental results show that the proposed model without initial contour can segment inhomogeneous images stably and effectively,which thereby avoiding the influence of contour initialization on segmentation results. Besides, the proposed model works better on noise images comparing with two relevant segmentation models.

Key Words: Image segmentation; partial differential equation; adaptive weight; local neighborhood;constant initialization.

1 Introduction

Image segmentation is a basic but significant problem in the field of computer vision and image processing,the goal of which is to separate the target from background. There are various kinds of images in the world,and they are always susceptible to some contaminations,such as inhomogeneity,noise,blur and so on,which make accurate segmentation still a challenging task. In this paper,we focus on images with intensity inhomogeneity.

Intensity inhomogeneity always occurs in medical images,especially in X-ray radiography,computed tomography(CT)and magnetic resonance(MR)images[1-3]. It often appears as intensity variation across an image, which results from technical limitations or artifacts introduced by the object being imaged. Popular global region-based models,such as Chan-Vese model[4],generally assume that image intensities are statistically homogeneous,i.e.,roughly a constant[2]. Thus,these piecewise constant(PC)models have difficulty to deal with intensity inhomogeneity. Segmentation of such inhomogeneous medical images usually requires intensity inhomogeneity correction as a preprocessing step[5].

In order to handle directly intensity inhomogeneity,some more sophisticated models than PC model are developed. Vese and Chan[6]and Tsai et al. [7]independently proposed two similar models,which are called piecewise smooth(PS)model.These methods can segment more general images,but they are computationally inefficient. To this issue,some local region-based models have been proposed[1, 2, 8-12]. For example, Li et al.[2]proposed a region-scalable fitting(RSF)active contour model(originally termed as local binary fitting(LBF)model[1]),which employs local region information and can deal with intensity inhomogeneity. However, the RSF model is computationally expensive since four convolutions need to be computed at each iteration for implementation. For more efficient segmentation,Zhang et al. [8]presented a local image fitting(LIF)model based on a constraint of the differences between the fitting image and the original image.Both models are only based on region mean information without considering region variance and thereby may lead to inaccurate segmentation. In[9],Zhang et al. exploited the local image region statistics and obtained a local statistical active contour(LSAC)model.The LSAC model can be directly applied to simultaneous segmentation and bias correction for 3T and 7T magnetic resonance images. For accurate medical image segmentation,an adaptive local-fitting-based active contour model is proposed in[10].

Local region-based models,however,to some extent,are sensitive to contour initialization. Dissimilar initial contours may lead to dissimilar segmentation results, even false results. This means they may be up against the problems of how and where to define initial contours. To allow more flexibility for contour initialization, some active contour models considering both global and local image information come into being.Wang et al. [3] presented a local and global intensity fitting (LGIF) model in a variational level set formulation. The LGIF energy is a linear combination of a local binary fitting(LBF)energy[1]and a global intensity fitting(GIF)energy[4],which has the form ofELGIF=(1?ω)ELBF+ωEGIF. The LGIF model shares the advantages of both the LBF model and the CV model. Thus, it is more robust to contour initialization, while handles intensity inhomogeneity efficiently. In fact,embedding global information into local region-based models can improve their robustness to contour initialization. Nevertheless,for such models [13-15], the parameter balancing the local and global information,such asωin the LGIF, must be adjusted in reference to initialization and noise level,which limits their practical applications[16].

In this study, we introduce a novel segmentation model for images with intensity inhomogeneity in the form of partial differential equation. First,we construct a global intensity fitting term and a local intensity fitting term. Then,a tradeoff function is defined to adjust adaptively the weight between two fitting terms. Finally, a weighted regularization term related to local entropy is used to smooth the evolution curve. Meanwhile,a distance regularization term is considered for stable level set evolution. The proposed model can segment images with intensity inhomogeneity without tuning the weight parameter. Furthermore,the level set function in the proposed model can be initialized as a constant function,which avoids troubles from the need of initial contour. Besides,the proposed model works better on noise images comparing with two relevant segmentation models.

The remainder of this paper is organized as follows. In Section 2, we briefly review the CV model and the RSF model, which are closely related to our model. Section 3 introduces and analyzes the proposed segmentation model. The implementation and experimental results are given in Section 4. This paper is summarized in Section 5.

2 Related works

2.1 CV model

The CV model[4]is a typical global region-based model,which approximates the original image through a piecewise constant function. Given a gray imageI:Ω?R2→R,letCbe a closed contour in the image domain Ω. The fitting energy of CV model is formulated as follows:

whereλ1andλ2are positive parameters,φis the level set function which divides Ω into two non-overlapping regions:Ωinside(C)={x:φ≥0}and Ωoutside(C)={x:φ<0}. By calculus of variations,c1andc2can be obtained by the following formulation:

Obviously,c1andc2are two constants that approximate the average intensities in Ωinside(C)and Ωoutside(C). For homogenous images,c1andc2always are optimal constants that minimize the global fitting energy(2.1). However,if the intensities inside or outside the contourC={x:φ(x)=0}are not homogeneous,such constants may be far away from original image data. As a result, the CV model has difficulty to segment images with intensity inhomogeneity.

2.2 RSF model

In order to handle intensity inhomogeneity,Li et al. [2]employed local image information and proposed the well-known RSF model.They introduced a kernel function and defined the following energy functional:

whereKσis the Gaussian kernel function with standard deviationσ, andfi(x)(i=1,2)are two smooth functions obtained by minimizing the energy functional (2.3) for fixed functionφ,which are expressed as

whereU(x,r)={y:‖x?y‖2≤r}andris the radius ofU(x,r).

The fitting functionfi(x)(i=1,2) make use of neighborhood information to fit the overall image gradually and then the RSF model can segment images with intensity inhomogeneity well. However,the RSF model may get stuck in local minima if the contour is not initialized appropriately.This makes the RSF model sensitive to contour initialization[16]. There is no simple answer that applies generally to date. A possible reason is that the fitting energy in(2.3)only takes the local image information into account.

3 The proposed model

3.1 Adaptive fitting energy

To address intensity inhomogeneity as well as allow more flexibility for contour initialization, in this study, we consider the following fitting energy with an adaptive weight functionω(x):

The adaptive weight functionω(x)plays a key role in our model;it can assign weights to two fitting terms adaptively according to the image feature. In the following, we describe this function in detail.

For an image pixel x∈Ω,we denote its 4-neighborhood,diagonal neighborhood and 8-neighborhood regions asU4(x,r),UD(x,r) andU8(x,r), respectively. Image intensity range of the whole image region and three neighborhood regions are represented asRΩ,RU4,RUDandRU8,which is computed by:

in which Max(:) and Min(:) are the maximum and the minimum intensity value in the given region. Clearly,RU4,RUD,RU8∈[0,255]andRΩ∈(0,255]for gray images if it is not a constant image.

Follow the definition of image local contrast ratio in [17], we define a novel local neighborhood contrast based onU4,UDandU8:

wherer1(x)is the local neighborhood contrast ofU8,whiler2(x)is the local neighborhood contrast related toU4andUD. We can see thatr1(x),r2(x)∈[0,1]and thenc(x)∈[0,1]. In Fig.1, we show the values ofr1(x),r2(x) andc(x) for an image with intensity inhomogeneity.

It can be seen from Fig.1 thatr1(x) has larger values at the object boundary and smaller values in other areas. Thus,r1(x) can detect object boundary accurately.r2(x)quantifies precisely the fluctuate of intensity value in the whole image,and reduced the influence of abnormal points, such as noise. The proposed local neighborhood contrastc(x)can not only indicate object boundary but also has some anti-noise property,which helps to improve the adaptability and accuracy of adaptive weight function.

To suppress the influence of regions with high contrast on the weight,we employ the nonlinear functiony=tanh(k·x)to adjust the value ofc(x). In this study,we choosek=3 and obtainθ(x)=tanh(3·c(x)),θ(x)∈[0,1). Inspired by[17],the adaptive weight functionω(x)is defined as:

whereb(x)=1?tanh(2·θ(x)).

3.2 Weighted regularization term

In order to guarantee the smoothness of evolution curveC, it is necessary to consider a regularization term for the proposed model. Generally, we can smoothCby penalizing its length[1-4,14,18,20]

or its weighted length[19,21]

whereg(|▽I|)=1/(1+|▽I|2)is an edge indicator function. In order to reduce the impact of noise,|▽I|is often replaced by|▽Iσ|=|▽(Kσ?I)|. In the smoothed versionIσ,noise is really suppressed,but the object edge is also blurred. This may result in inaccurate edge indication.

The value ofg(E) decreases when the local entropyEincreasing, andg(E)→0 whenE→∞. The weighted regularization term is finally formulated as

For a point x∈Ω and its neighborhood Ωkwith sizeMk×Nk,the local entropy of the neighborhood Ωkis defined by[22]:

whereLis the maximal gray level of imageI,pi=ni/(Mk×Nk) is the probability of grayscaleiappears in the local region Ωk,niis the number of pixels with grayscaleiin Ωk. Local entropy reflects the variation of image intensity in a neighborhood. The local entropy value is large if the intensity in the neighborhood is homogeneous, while the value is small for a heterogeneous or noisy region.

3.3 Total evolution equation for level set function

As in typical level set methods,we introduce a distance regularization term[23]into our model for stable level set evolution:

With the above defined fitting energy term(3.2),the weightedL2regularization term(3.10)and the distance regularization term(3.12), our model is finally formulated as the following partial differential equation:

with the initial and Neumann boundary conditions:

In most of the level set models, the level set function is initialized to be a signed distance function [4, 6] or a binary function [1, 2, 23]. The proposed model allows for more flexible contour initialization. Especially,the initial level set function can be chosen as a constant function:

whereρis a positive constant. Such constant initialization of level set function not only facilitates to operate in practice, but also effectively avoids the sensitivity to initial contour. We do not need to consider the problems such as where and how to initialize the contour.

4 Implementation and experimental results

4.1 Numerical implementation

To compute the evolution equation (3.13) numerically, we choose slightly regularized version of the Heaviside functionH(·)and Dirac functionδ(·)as follows[4]:

withε=1.0.

The evolution equation(3.13)can be implemented by simple finite difference method.Lethstand for the spatial step of the discrete grid,and △tbe the time step. The level set functionφ(x,t)=φ(x,y,t) can be discretized asφni,j=φ(ih,jh,n△t). The temporal partial derivative is approximated by forward difference, while two spatial partial derivatives?φ/?xand?φ/?yare approximated by central difference. The Laplace operator △φis calculated by basic five-point difference scheme. Then, the approximation of (3.13) can be simply written as:

in whichA(φn+p) andA(φn) represent the area enclosed byC={x:φ(x)=0} with an interval ofpiterations,ξis a very small positive real number,and we takep=10,ξ=0.01 in the following experiments.

The procedures of the algorithm for solving problem(3.13)can be summarized in Algorithm 4.1.

Algorithm 4.1:1. Initialize φ0=2,the maximum iteration number N=500,and set n=0;2. Calculate ω(x)using(3.6);3. Compute c1(φni,j)and c2(φni,j)by(2.2),and compute f1(φni,j)and f2(φni,j)by(2.4);4. Update φni,j according to(4.3)and obtain φn+p i,j ;5. Compute A(φn+p)and A(φn),then check whether the termination criterion(4.4)is satisfied or n+p≥N. If not,let n=n+1 and return to step 3.

4.2 Experimental results

In this section, the proposed model is applied on some synthetic images as well as real images, and is compared with three relevant models: RSF [2], LGIF [3] and LSAC [9].We show experimentally that the proposed scheme can segment images with intensity inhomogeneity effectively and be more robustness to noise without tuning the weight between two fitting terms.In addition,in our model,the level set function can be simply initialized to a constant function, which avoids inherent problems caused by contour initialization.

All experiments presented in this section were performed on an Intel Core i7-2600,CPU 3.40GHz processor. The models with algorithms were implemented in MATLAB(R2016a)software.Unless otherwise specified,we use the following parameters:△t=0.1,λ1=λ2=1.σandλare adjusted according to different kinds of images. In general,if the image has more details and its intensity varies very large, the parameterσshould take a smaller value. The parameterλis chosen according to noise level, i.e., larger values for images with heavy noise. Their default values are setσ=3,λ=0.0005×2552. The parameters for RSF,LGIF and LSAC models are taken from the associated papers[2, 3,9]. The values of weightωfor the LGIF model and the Gaussian kernelζfor the LSAC model are given in the experiment.

Firstly, in Fig.2, we test RSF, LGIF, LSAC and our models on four typical images with intensity inhomogeneity. The first row of Fig.2 displays the original images as well as initial contours for RSF, LGIF and LSAC models. The segmentation results of RSF, LGIF and LSAC models with initial contours are shown in the second, third and fourth rows, respectively, while the segmentation results of our model without initial contour are presented in the last row. In this experiment, we tuned parameters for optimal results of all models and finally chooseν=0.008×2552for the RSF model,ν=0.008×2552,ω=0.03,0,0.1,0.1 for the LGIF model,ζ=3,6,6,1 for the LSAC model andλ=0.0006×2552,0.00035×2552,0.00052×2552,0.0008×2552for our model. For Fig.2(a)-(c),four models all obtain accurate segmentation results. For Fig.2(d), the LSAC model fails to segment it(Fig.2(p))and other three models can still obtain satisfying segmentation results.

The corresponding CPU times of four models for segmenting images in Fig.2 are listed in Table 1. The odd rows of Table 1 give the CPU times for the whole segmentation,while the even rows display the average CPU times for single iteration. We observe that,without initial contour,our model can not only handle intensity inhomogeneity effectively, but also has less segmentation time among four models. Note that, compared with other three models, the proposed model needs more times in average for single iteration due to the calculation of the weight of all image pixels and the termination criterion(3.18). However,the proposed model greatly reduces iteration numbers to obtain segmentation results thanks to the adaptive weight function.

Secondly, we test the sensitivity of the RSF, LGIF and our models to the position of initial contour in Fig.3. Due to the space limitation,only an X-ray vessel image is chosen as the representative test image,which is a typical image with intensity inhomogeneity,as shown in Fig.3(a). We choose five different initial contours with the same size but different positions, as shown in the first row of Fig.3. From left to right in Fig.3, the range of weightωin LGIF is [0.048,0.074], [0.001,0.066], [0.036,0.075], [0.082,0.177] and[0.039,0.089], respectively. For the proposed model,we chooseλ=0.0007×2552andσ=9. The segmentation result of our model without initial contour is listed in Fig.3(b),while the results of the RSF and LGIF models with these initial contours are listed in the second and third rows. From Fig.3, we can see that the RSF model is sensitive to the position of initial contour and can only segment the vessel image for the last initial position(Fig.3(l)). The LGIF model can reduce its sensitivity by manually adjusting the weightω,but there are still some unsatisfactory results,as shown in Fig.3(o)and(p).By contrast,the proposed model accurately extracts the object without initial contour.

Figure 1: Example of the local neighborhood contrast. (a): test image. (b): 2-D plot of r1(x). (c): 2-D plot of r2(x). (d): 3-D plot of c(x). (e): 3-D plot of r1(x). (f): 3-D plot of r2(x).

Image ID RSF model LGIF model LSAC model proposed model(a)10.3763 s5.0485 s10.0986 s1.1800 s 0.0158 s0.0154 s0.0281 s0.0179 s(b)10.1945 s11.2114 s7.9786 s2.2770 s 0.0244 s0.0269 s0.1308 s0.0278 s(c)3.7305 s3.1464 s5.8428 s0.6526 s 0.0178 s0.0178 s0.0860 s0.0186 s(d)14.1013 s8.8999 s33.9217 s5.3447 s 0.0644 s0.0685 s0.0424 s0.0722 s

Thirdly,we compare the robustness of the RSF,LGIF,LSAC and our models to noise.Three noise images together with initial contours are shown in the first column of Fig.4.In this experiment, we also adjusted parameters for optimal results of all models and finally chooseν=0.008×2552,0.01×2552,0.003×2552in the RSF model,ω=0.42,0.5,0.01 in the LGIF model,andν=0.05×2552,0.02×2552,0.001×2552,ζ=6 in the LSAC model.For our model,we setλ=0.0015×2552,0.0016×2552,0.0007×2552.

From the segmentation results shown in the second and third columns in Fig.4, we observe that the RSF model and the LGIF model successfully extract the object for the first image (Fig.4(b), (c)), but fail to segment other two images (Fig.4(g), (l), (h), (m)).The segmentation results of the LSAC model are shown in the fourth column, in which we can see that it is easily susceptible to noise. The proposed model provides accurate segmentation for all images without initial contour,as shown in the last column in Fig.4.

Finally,in Fig.5,we apply the proposed model on four images with different modalities. The original images are listed in the upper row and the segmentation results of the proposed model are displayed in the bottom row of Fig.5. Fig.5(a)is a real plane image with uneven illumination and Fig.5(b)is a real chip image with complicated background.From Fig.5(e)and(f),we can see that the proposed model has successfully extracted the edges of objects from both images. The third image shown in Fig.5(c) is a multiphase synthetic image with weak and blurry edges. The proposed model can detect all objects simultaneously and accurately with only one level set function. The last image listed in Fig.5(d)is a very different one: range image of saddle with step or roof edges,which is taken from[24]. This range image has complicated image intensity distribution,in which the end boundaries of the saddle are step edges,and the side boundaries are roof edges.Because the image intensity varies throughout the object corresponding to the distance from the viewer,it is difficult for the constant approximation to detect and segment this kind of object [24]. As seen from Fig.5(h), the proposed model detects both step edges and roof edges exactly.

Figure 2: Segmentation for four typical images with intensity inhomogeneity. Row 1: initial contours. Row 2:results of the RSF model. Row 3: results of the LGIF model. Row 4: results of the LSAC model. Row 5:results of the proposed model.

5 Conclusion

In this study,we define a novel neighborhood contrast based on the neighborhood of an image point,and derive an adaptive weight function to adjust the weightiness between global and local fitting terms. Besides,a weighted regularization term related to image local entropy is introduced to smooth the evolution curve. Combined with a distance regularization term, our model is finally formulated in the form of partial differential equation. The proposed model is not only effective in segmenting images with intensity inhomogeneity, but also robust to noise. Especially, the proposed model does not need to set initial contour,which essentially avoids some problems associated with the initial contour. However,when the model is applied to large images,it is time consuming due to the calculation of local neighborhood contrast. In the future research,we will focus on the fast algorithm as well as its combination with new segmentation methods.

Figure 3: Comparison of three models for an X-ray vessel image in terms of robustness to initial contour position.Row 1: original image and result of the proposed model without initial contour. Row 2: initial contours. Row 3: results of the RSF model. Row 4: results of the LGIF model.

Figure 4: Applications of the RSF, LGIF, LSAC and our models on three noise images. Column 1: initial contours. Column 2: results of the RSF model. Column 3: results of the LGIF model. Column 4: results of the LSAC model. Column 5: results of the proposed model.

Figure 5: Segmentation of the proposed model on images with different modalities. Row 1: original images.Row 2: results of the proposed model.

Acknowledgments

The authors would like to thank anonymous reviewers. This work is partially supported by National Natural Science Foundation of China (No. 11901071 and 31971113), the Natural Science Foundation Project of CQ CSTC (No. cstc2019jcyj-msxmX0219), the Educational Commission Foundation of Chongqing of China (No. KJQN201800537), the China Postdoctoral Science Foundation (No. 2020M671685) and the Selective Funding for Postdoctoral Research Projects in Zhejiang Province(No. 514000-X81902).