International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 56

ISSN 2229-5518

Level Set Method for Image Segmentation

and Bias Field Estimation using Coefficient of

Variation with Local Statistical Information

Mr.Said Anwar Shah, Dr. Noor Badshah**,**

Abstract— Most of the image processing techniques use image regional information for image segmentation, image registration, etc. Region based methods for image segmentation with bias field estimation based on the statistical information of different region such as intensity mean, intensity distribution etc. These methods rely on the assumption of spatial invariance in the image domain that is not always valid for most images. We present a new variational model in this paper for segmenting important tissues in the image with the estimate of bias field. We use a term like coefficient of variation as a criterion function based on the local statistical information, and define energy functional in the level set framework. Coefficient of variation is a relative measure uses relative intensity information in the neighborhood for achieving the desired result. Minimization of this energy functional gives estimate of bias field and segment important tissue in the image domain. Empirical results reveal that our new proposed model is independent of initialization and accurate in contrast with the existing model.

Index Terms—Image segmentation, Intensity inhomogeneity, Coefficient of variation (CoV), Bias estimation, level sets, MRI, Functional minimization.

—————————— ——————————

spatial invariance; example is Chan and Vese model [1].

This model does not work for the images with

RI is a powerful imaging techniques in medical images used for the visualization of inner body

structure, while bias field is unwanted smooth signal that occur in most MRI due to intensity inhomogeneity in their magnetic field. This undesirable smooth signal modifies high frequency information of the image, which create many problems in image processing algorithms such as image segmentation, image registration, image classification etc. In the presence of bias field it is difficult to segment important tissue in the image domain due to intensity overlapping in different region. Most of the image segmentation algorithms [1],[2] rely on the assumption of spatial invariance, which are not applicable to images with bias field. Level set methods as a numerical techniques used for the evolution of contour by embedding the contour in one higher dimension is called level set function. This method is design for a problem to move the contour front forward in some direction if the

heterogeneous regions. Further Chan and Vese generalized his work [3] into multiphase level set formulation by using multiphase level set functions for different regions such models are called piecewise smooth (PS). These models do not base on the assumption of intensity homogeneity and work well even in the presence of intensity inhomogeneity. However, it is difficult to compute numerically.

2 BACKGROUND REVIEW

2.1 Mumford and Shah Functional Model

Among variational models Mumford and Shah model is the oldest model used for image segmentation. Let *v *and *I *are the functions define on Ω ⊂ R 2 .The aim of image segmentation is to find out the desired *v *from the given image *I *through a defined model. Mumford and Shah proposed the following variational model [6] gives disjoint

region in the image domain

level set function is positive otherwise backward to reach

the desired object boundary. Active contour methods commonly use two approaches; Edge based models [4],[ 5]

*E*(*v*, *C *) = λ1 ∫ | *v *− *I *|

Ω

dxdy +λ2

∫ | ∇*v *|2

Ω \*C*

*dxdy *+λ3 | *C *| ,

(1)

and Region based models [1],[2],[ 3] for classifying pixel

where λ1 > 0, λ2 > 0

and λ3 > 0 are parameters use for

value belonging to object or region. Edge-based models uses important feature in the image domain for segmentation such as edges that separate regions. These models even works in the presence of bias field, but it is hard to find edges of different regions in noisy images. Region-based models are superior to edge-based model, because it uses regional information for image segmentation that take more pixel value than edge and will have more information for the movement of active contour. Region-based models rely on the assumption of

the adjustment in the model, and | *C *| gives the boundary

length. Eq.1 has three terms. First term is data fitting term, second term use to regularize the contour, and third is the smoothing term. The minimization of this energy functional through partial differential equation has some problems. The main problem is discrete values on boundaries, which create discontinuities into energy functional. To surmount this problem Chan-Vese [1] used level set methods and solved Mumford and Shah Segmentation Problem.

IJSER © 2014 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 57

ISSN 2229-5518

3 Framework for Image Segmentation with

intensities values in O y ∩ Ω i

form a cluster with center

mi =b(y)

ci , may be treated are sample values drawn

Image processing techniques for image segmentation uses the image model for segmenting important region in the image domain. It assumed to have two important components, the bias field and the true image. We use

multiplicative model for modeling bias field, an observed

from Gaussian distribution with mean mi .

In the neighborhood O y we use local intensity clustering property to classify intensity values into N cluster with

image *I *can be expressed in term of true image *q*, bias field

centers

mi =b(y) ci

i=1, 2, 3... N. K-mean clustering

b, and additive noise n as

I = bq + n,

(2)

algorithm minimizes the clustering criterion function

through iterative process [7]. In the proposed model, we

where *b *is the bias field, *q *is the true image presumed to be constant in different region and *n *assumed to be zero mean additive Gaussian. Our basic objective is to extract bias field components from the given image and use them to correct intensity non uniformity. We assume that bias field is smooth and slowly varying throughout the domain

use a modified form of K-Mean clustering criterion

function to classify the intensity values in the

neighborhood O y into disjoint regions. For the intensities *I*(x) the continuous version of the clustering criterion function can be written as

N

can be approximated by constant values in the

neighbourhood. The true image *q *will have constant

R y = ∑

i =1

∫ | *I *( *x*) − *m*i |

Oy

vi dx

(5)

values *c*1 , *c*2 ,...*c*n in disjoint region Ω1 , Ω 2 ,...Ω n reflect

where vi is the membership function of the region Ω i *,*

the physical properties in these regions.

vi =1 for x ∈ Ω i and vi =0 for x ∉ Ω i , mi

is the ith

Region based methods use different region descriptor of

cluster center. The above energy equation can be rewritten as

N

the intensities, such as intensity mean, intensity distribution, etc to segment important tissue in the image

R y = ∑

i =1

∫ | *I *( *x*) − *m*i | *dx*

Ωi ∩*O*y

(6)

domain. In the presence of bias field it is to define a region

The clustering criterion function in (6) is defined in term of

descriptor that gives information about different regions in

the image domain, because of overlapping between the intensity distributions in Ω1 , Ω 2 ,...Ω n regions. So, it is

centre approximation by intensities in O y as

N

mi =b(y) ci

for classifying the

not possible to use intensity information for region

segmentation. However, in the formulation of our proposed model we use the local intensity clustering

R y = ∑

i =1

∫ | *I *( *x*) − *b*( *y*)*c*i

Ωi ∩*O*y

|2 *dx*

(7)

property for achieving the desired result. According to this property at every point y belong to Ω , we define a circular neighborhood with radius *r *is given by

We introduce a nonnegative function window function

K(y-x) is defined as K(y-x)=0 for x ∉ O y . The energy function in term of center approximation and window

{*O *y

= *x *:| *x *− *y *|< *r*}.

(3)

function can be rewritten as

N

Partition of the full domain Ω also partition of the small

N

Fy = ∑

∫ K ( y − x) | I ( x) − b( y)ci | dx

(8)

circular region O y

i.e.

O y ∩

Ω form partition O y .

i =1 Ωi

i =1

The values of the smooth function b(x) for x ∈

O y are

Minimization of the energy functional in Eq.8 results

image segmentation with bias field estimation. It does not

close to b(y) i.e. b(x) ≈ b(y) for x ∈

O y . Therefore the

work well in the images with more intensity Non- uniformity. Therefore, in the proposed model we present a

intensity value b(x) q(x) in subregion O y ∩ Ω i are close to

simple modification of the energy functional in Eq.8

b(y) ci

for x ∈ *O *y ∩ Ω i

i.e. *b(x)q(x) ≈ b(y) *ci

for x ∈

dividing the criterion function by the square of the

O y ∩ Ω i

where *q*(x) represent true value and function of

approximated cluster mean b(y) ci

to exploit the relative

x and *b(y) *ci are the estimated constant values in different

information of the cluster. The modified form of energy in

Eq.8 can be defined as

regions. Thus, the observed image in (2) can be rewritten N

| *I *( *x*) − *b*( *y*)*c*

|2

as ∑ ∫

i

I = b( y)ci + n(x),

(4)

Fy =

i =1

K ( y − x)

Ωi

(*b*( *y*)*c*i )

*dx*

(9)

where n(x) Gaussian noise with zero mean and the

Our proposed method uses the modified form of the K-mean

clustering criterion function as tool in the energy functional F

y

IJSER © 2014 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 58

ISSN 2229-5518

to classify the intensities in the neighbourhood O y , which

N

Ω 2 *, *a single level set function φ is used in the level set

formulation. The membership functions such as

partition O y such as

O y ∩

. Smaller the value of the

i

*D*1 (φ ) = *H *(φ ) ) and

*D*2 (φ ) = 1 − *H *(φ )

are used to

i =1

define disjoint regions

Ω1 and

Ω 2 . Where H is the

energy functional

Fy result good classification. Good

Heaviside function [8],[9]. Thus, the energy functional in

(10) in term of level set formulation can be written as

classification of the intensity values in the entire domain Ω can 2

| *I *( *x*) − *b*( *y*)*c*

|2

be obtained by minimizing

Fy for all value of y in Ω . Thus we

i

i

(*b*( *y*)*c *) 2

*D*i (φ ( *x*))*dx **dy*

minimize

Fy for all value of y in Ω , which can be obtained

i =1 Ω

(11)

minimizing ∫ Fy dy for all y in Ω such as F = ∫ Fy dy

By exchanging the order of integration we have,

N | *I *( *x*) − *b*( *y*)*c*

|2

| *I *( *x*) − *b*( *y*)*c*i |

i

F =

K ( y − x)

*dy **D *(φ ( *x*))*dx*

*F *= ∫ ∑ ∫ *K *( *y *− *x*)

(*b*( *y*)*c *) 2

*dx **dy*

(10)

∫ ∑ ∫

(*b*( *y*)*c *) 2 i

i =1 Ωi

i

i =1 Ωi

i

In Eq.10 the integration is over the entire domain, therefore, we omit Ω i in the subscript. Minimization of the energy functional in (10) with respect to the regions,

(12)

The above energy *F *depends on the level set function φ ,

the constants c1 , c2 ,...cn , and the bias field b can be

Ω1 , Ω 2 ,...Ω n

, bias field *b*, and constants

c1 , c2 ,...cn

rewritten as

segment important tissues in the images with estimate of 2

bias field. In the proposed method we use truncated

*F *(φ , **c**, *b*) = ∫ ∑ *a*i ( *x*)*D*i (φ ( *x*))*dx *

(13)

Gaussian distribution as a kernel function defined by

i =1

where ai is given by

1

K (v) =

− 1 (|*v*|/ σ )2

e 2

for

all

| *v *|< *r*

| *I *( *x*) − *b*( *y*)*c*

i

|2

0

otherwise

ai ( x) = ∫ K ( y − x)

(*b*( *y*)*c *) 2

*dy *. (14)

where ϭ is the scale parameter of the Gaussian function, β

The simple expression for the function ai

is given by

is the normalization constant and *r *is radius of the

2 1 *I *2 1

neighbourhood

O y . We select *r *for the neighbourhood

ai ( x) = I k −

ci

I ( * K ) + ( * K ) , (15)

b c 2 b 2

O y on the basis of the degree of intensity inhomogeneity.

The bias field vary faster in more intensity inhomogeneous regions; therefore, we take small value for the radius *r *in the neighbourhood O y .

where * is the convolution operator, and *I *k = ∫ *K *( *y *− *x*) is one in the image domain Ω except near the boundary. The data fitting term in (13) can be expressed in a variational level set framework is given by

4. Energy Minimization using level set

*R*(φ , *b*, **c**) = *F *(φ , *b*, **c**) + λ1 *R*q (φ ) + λ2 *L*(φ )

(16)

where

*R*(φ ) is the distance regularize have been used in

We use level set method that was devised by Osher [7]. In the propose method we use level set functions takes

different variational method [8] for the smooth curve

evolution is given by

positive and negative signs to represent different region in

*R*q (φ ) = ∫ *q*(| ∇φ | *dx*

(17)

the image domain. It is difficult to minimize the above

energy in the present form; therefore, we use level set

formulation and well known variational method for the

and *L*(φ ) is the arc length of the contour can be defined

as

minimization of energy in (10). Suppose

φ : Ω → R be a

*L*(φ ) = ∫ | (∇*H *(φ ) | *dx*

(18)

level set function, its signs partition the domain into two

The minimization of the energy in (16) results image

segmentation with bias field estimation and this can be

disjoint regions such as Ω1 = {x: φ (x)>0} and

Ω 2 = {x:

obtained by iterative process. For minimization, we

φ (x) <0}. It is called two phase level set formulation. To

partition the image domain Ω into more than two

partially differentiate the energy functional *R*(φ , *b*, **c**) in every iteration with respect to φ ,b and c.

disjoint regions Ω1 , Ω 2 ,...Ω n

, we use two or more than

two level set functions are called multiphase level set formulation.

To partition the image domain into two regions Ω1 and

IJSER © 2014 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 59

ISSN 2229-5518

The energy functional minimization in (16) with respect to

φ can be obtained by using the Gateaux derivative ∂R of

∂φ

∧

iii) Update *b *to be b for fixed value of φ and **c.**

∧

iv) Update * c *to be c for fixed value of φ and

We implement our model by using finite difference

the functional *R*(φ , *b*, **c**) , and solve the following gradient

scheme with time step δt=0.1, parameters

λ1 =1 and

flow equation for the steady state solution.

λ2 =0.001*255^2. We set the scale parameter ϭ of the

∂φ = − ∂*R*

∂*t *∂φ

(19)

Gaussian distribution is 4 and keeping small size of the neighbourhood O y for the images having more intensity

The expression for the Gateaux derivative can be obtained by using calculus of variation [9], hence the gradient flow equation is given by

∂φ

inhomogeneity. Fig.1 shows the result for MRI images. The curve evolution processes has given in the images by showing the initial contours on the images (in the first

column), estimated bias field (in the second column) and

= −δ (φ )(*a*1 − *a*2 ) + λ1 *div*(*d*

∂*t *q

(| ∇φ |)∇φ )

∇φ

(20)

bias corrected images (in the third column).

+ λ2δ (φ )*div*(

)

| ∇φ |

Finite difference method can be used to implement the smooth evolution of level set function in (20) and update

the constant value

(*c*1 , *c*2 )

in **c **and *b *through

minimization of energy functional with respect to **c **and b respectively.

Minimize the energy functional *R*(φ , *b*, **c**) with respect to

**c **keeping φ and b are constant gives the optimal value of

c= (*c*1 , *c*2 ) as

1

∧ ∫ 2

* K )Ivi dy

c = __ b __

( 1 * *K *)*v dy*

i = 1,2, with vi

= *D*i (φ ( *x*) (21)

There is an obvious intensity inhomogeneity in the above

∫ b i

4.4 Energy Functional Minimization With Respect to b:

images. Our method exploits relative intensity information of the images to obtain the estimate unwanted bias field. We correct intensity inhomogeneity by using the estimated value of bias field b and bias corrected image can be

Minimize the energy functional *R*(φ , *b*, **c**)

with respect

obtained dividing the image I by estimated bias field b. To

to b keeping φ and c are constant gives the optimal value of b as.

know about the quality of our proposed method, we test

our method on medical images with clear intensity

∧ (W

b =

( 2)

I * K ) ,

(22)

inhomogeneity. In Fig.2 the initial contour is shown on the

original image in the first column and their corresponding

W (1) * K

Where,

segmentation result, bias field estimation, and bias

corrected image are shown in second, third, and fourth

2 (1)

2 ( 2)

column, respectively.

W = ∑

i =1 i

vi ,

and W

= ∑ 2 i

i =1 ci

We use the following algorithm for the minimization of energy functional in (16) as given below:

i) Initialization of φ ,b and **c.**

∧

ii) Update φ to be φ for fixed value of *b *and **c.**

IJSER © 2014 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 60

ISSN 2229-5518

It can be seen in Fig.2 the intensities within different region are heterogeneous in the original images and become homogeneous in the bias corrected images. We can judge the result and quality of our proposed model through histogram analysis by observing the histogram of the uncorrected and bias corrected images. The histogram of the original images not showing well- defined and well- separated peak due to the presence of bias field in their intensity distribution. The uncorrected image with initial contour, estimated bias field, bias corrected image are shown in first, second, and third column of row1, histogram of the uncorrected and bias corrected images are shown in the first, second column of row2 in Fig.3.

It can be seen in Fig.3 the histogram of the bias corrected images have well-separated peaks shows intensity homogeneity of different region, while histogram of the uncorrected image are not having well-separated peaks due the presence of bias field.

Since the Chunming Li model perform well as compared to other active contour models in the presence of bias field, we compare the result of our model with the Chunming Li (CLIC) model in [10]. For this purpose, we are using histogram analysis, because histogram of the bias corrected images will have well- separated peaks showing different homogenous regions in the image domain. Histogram in the second row of Fig.4 is the result of our proposed model having well-defined and well-separated peaks showing different homogenous regions in the bias corrected image, while histogram in the fourth row of Fig.4 is the result of Chunming Li model in [10] do not have well-separated peaks due to the presence bias field.

5. CONCLUSION:

Bias field is a low frequency signal that modifies the intensity value of the images; therefore it is necessary to reduce its value before any image processing techniques. In this paper, we have presented a novel approach for image segmentation with bias field estimation work well

IJSER © 2014 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 61

ISSN 2229-5518

in the presence of intensity inhomogeneities. We have used local intensity clustering property and define energy functional like coefficient of variation in term of level set function, which uses relative intensity information in small neighbourhood. Minimization of this energy functional result jointly image segmentation with bias field estimation. Empirical results show that our proposed method produce best result in term of efficiency, accuracy and independent of initialization. .

[1] T. F. Chan and L. A. Vese, Active contour without edges.

IEEE Trans. Image process, 10(2): 266-–277,2001.

[2] R. Ronfard, Region-based strategies for active contour model. J.

Comput. Vis, 13(2):229—251,1994.

[3] L. Vese, T. Chan, A multiphase level set framework for image segmentation using Mumford and Shah Model. J. Comput. Vis,

50(3):271--293, 2002.

[4] V. Caselles, R. Kimmel, and G. Sapiro, Geodesic active contours, J. Comput. Vis, 22(1): 61–-79, 1997.

[5] R. Malladi, J. A. Sethian, and B. C.Vemuri, Shape modeling with front propagation: A level set approach, IEEETrans. Pattern Rec. 17(2): 158–175,1995.

[6] D. Mumford and J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems. Commun.Pure Appl. Math, 42(5): 577–685, 1989.

[7] S. Osher and J. A. Sethian, 'Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79, 12–49 ,1988

[8] C. Li, C. Xu, C. Gui, and M. D. Fox, 'Distance regularized level set evolution and its application to

Image segmentation, IEEE Trans. ImageProcess, 19(12): pp.

3243–3254, 2010.

[9] K. Zhang, L. Zhang, H. Song and W. Zhou, Active Contours with selective local or global segmentation: A new formulation and level set method, Image and vision Computing, 28 : 668--

676, 2010.

[10] C. Li, R.Huang, Z.Ding,C. Gatenby,D. N. Metaxas, and J. C.

Gore," A Level Set Method for Image segmentation in the presence of Intensity Inhomogeneities With application to MRI" IEEE Trans. ImageProcess, 20(7) :2007-2016, 2011.

• *Mr. Said Anwar Shah is currently pursuing(MS) degree program in*

Mathematics at University of Engineering and Technology Peshawar

,Pakistan., PH-091-9216502, E-mail: anwarshahuet@yahoo.com

• *Dr. Noor Badshah Assistant Professor department of Basic Science and*

Islamiat UET Peshawar, Pakistan, Ph-091-9216502,E- mail: noor2knoor@gmail.com..

IJSER © 2014 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014

62

ISSN 2229-5518

I£ER lb) 2014