[Journal of mathematics at home University] No-edge motion Contour Model in 054th Image Segmentation

$ \ BF abstract $: This article gives a detailed description of section 4.4 of "methods of partial differential equations in image processing" Edited by Wang dakai.

$ \ BF keywords $: Image Segmentation; Active Contour Model; matlab programming

 

1. Model Creation

In an image, the difference between the object and the background is sometimes obvious difference in average gray scale. because such images do not have obvious edges ($ \ sev {\ n I} $), nor do they have obvious textures (texture and gray scale changes have a certain pattern, therefore, the Geodesic Active Contour, GAC, or snake model cannot achieve image segmentation.

For this reason, T. chan and L. vese proposed the following energy functional: $ \ bee \ label {1: e} e (C_1, C_2, c) \ equiv \ MU \ oint_c \ rd s + \ lambda_1 \ iint _ {\ omega_1} (I-c_1) ^ 2 \ RD x \ RD y + \ lambda_2 \ iint _ {\ omega_2} (I-c_2) ^ 2 \ RD x \ RD y, \ EEE $ where

(1) $ \ Mu> 0 $, $ \ lambda_1 $, $ \ lambda_2> 0 $ is the weight of each score;

(2) The curve $ C $ is simple, and the image area $ \ Omega $ is divided into two parts, respectively, as $ \ omega_1 $, $ \ omega_2 $ ;\

(3) $ C_1 $, $ C_2 $ is a two-constant, table shard constant image.

This model is characterized

(1) On the one hand, shortest the boundary curve;

(2) On the other hand, make the background $ (\ omega_2) $ different from the object $ (\ omega_1) $ (the last two items in \ eqref {1: e ). this is called an active contour without edge, C-V, or Gar model.

Introduce the $ \ Delta $ function and write $ e $ as $ \ bee \ label {1: EU} \ Bea E (C_1, C_2, U) & =\ MU \ iint _ \ Omega \ delta (u) \ sev {\ n u} \ RD x \ rd y \ & \ quad + \ lambda_1 \ iint _ {\ omega_1} (I-c_1) ^ 2 \ SEZ {1-H (u)} \ RD x \ RD y + \ lambda_2 \ iint _ {\ omega_2} (I-c_2) ^ 2 H (u) \ RD x \ RD y, \ EEA \ EEE $ where

(1) $ \ delta (\ cdot) $ is a Dirac $ \ Delta $ function with smooth approaching elements $ \ Bex \ Delta _ \ ve (X) = \ frac {1} {\ PI} \ cdot \ frac {\ ve} {x ^ 2 + \ ve ^ 2} \ quad \ sex {0 <\ ve \ LL 1 }; \ EEx $

(2) $ \ sev {\ n u} $ is the transformation of the Jacob;

(3) $ H (\ cdot) $ is the Heaviside function: $ \ Bex H (z) =\left \{\ BA {ll} 1, & Z \ geq 0, \ 0, & Z <0, \ EA \ right. \ EEx $ it has smooth approaching yuan $ \ Bex H _ \ ve (z) =\ frac {1} {2} + \ frac {1} {\ PI} \ arctan \ frac {z} {\ ve }, \ EEx $ and $ \ Bex H' (z) = \ delta (z), \ quad \ mbox {In} \ mathcal {d} '(\ BBR) \ mbox {top }. \ EEx $

(4) under the $ U $ fixed condition of the function, by $ \ Bex \ iint _ {\ omega_1} (I-c_1) ^ 2 \ RD x \ RD y = C_1 ^ 2 \ iint _ {\ omega_1} \ RD x \ RD y-2c_1 \ iint _ {\ omega_1} I \ RD x \ RD Y + \ iint _ {\ omega_1} I ^ 2 \ RD x \ RD y, \ EEx $ \ Bex \ iint _ {\ omega_2} (I-c_2) ^ 2 \ RD x \ RD y = C_2 ^ 2 \ iint _ {\ omega_2} \ RD x \ RD y-2c_2 \ iint _ {\ omega_2} I \ RD x \ RD Y + \ iint _ {\ omega_2} I ^ 2 \ RD x \ RD y, \ EEx $ and the properties of quadratic functions are known only when $ \ bee \ label {1: c1c2} c_1 = \ frac {\ iint _ {\ omega_1} I \ RD x \ RD y} {\ iint _ {\ omega_1} \ RD x \ RD y }, \ quad C_2 =\frac {\ iint _ {\ omega_2} I \ RD x \ RD y }{\ iint _ {\ omega_2} \ RD x \ RD y} \ EEE $ $, $ e $ get minimum.

Now, write the variation (variation) of \ eqref {1: EU) $ \ Bex \ frac {\ Delta e} {\ Delta u} =-\ MU \ delta (u) \ Div \ sex {\ frac {\ n u} {\ sev {\ n u}-\ lambda_1 \ delta (u) (I-c_1) ^ 2 + \ lambda_2 \ delta (u) (I-c_2) ^ 2. \ EEx $

The gradient descent flow is $ \ bee \ label {1: GDF} \ frac {\ P u} {\ p t} =-\ frac {\ Delta e} {\ Delta u} = \ MU \ delta (u) \ Div \ sex {\ frac {\ n u} {\ sev {\ n u }}+ \ lambda_1 \ delta (u) (I-c_1) ^ 2-\ lambda_2 \ delta (u) (I-c_2) ^ 2. \ EEE $

2. Numerical Algorithms

For ease of calculation, discretization is required, and $ \ Delta _ \ ve $ is used to replace $ \ Delta $. in this way, the semi-implicit (semi-implicit) solution is $ \ bee \ label {2: scheme} U ^ {n + 1 }_{ IJ} = u ^ N _ {IJ} + \ Tau \ Delta _ \ ve (U ^ N _ {IJ }) \ SEZ {\ Mu Q (U ^ {n + 1 }_{ IJ}) + (I _ {IJ}-C ^ N_1) ^ 2-(I _ {IJ}-C ^ N_2) ^ 2}, \ EEE $ selected $ \ lambda_1 = \ lambda_2 = 1 $, $ q (U ^ {n + 1 }_{ IJ}) $ adopts the format of combining forward and backward difference: $ \ Bex Q (U _ {IJ }) & = d ^ {(+) }_x \ sex {g _ {IJ} d ^ {(-)} _ x u _ {IJ }}+ d ^ {(+)} _ Y \ sex {g _ {IJ} d ^ {(-)} _ Yu _ {IJ }}\ quad \ sex {g _ {IJ }=\ frac {1 }{\ sev {\ n u }_{ IJ }}}\\& = g _ {I, J + 1} d ^ {(-)} _ Xu _ {I, j + 1}-G _ {I, j} d ^ {(-)} _ Xu _ {IJ} + G _ {I + 1, J} d ^ {(-)} _ Yu _ {I + 1, j}-G _ {IJ} d ^ {(-)} _ Yu _ {IJ} \ & = g _ {I, J + 1} \ sex {u _ {I, J + 1}-U _ {IJ}-G _ {IJ} \ sex {u _ {IJ}-U _ {I, j-1 }}\& + G _ {I + 1, J} \ sex {u _ {I + 1, j}-U _ {IJ}-G _ {IJ} \ sex {u _ {IJ}-U _ {I-1, j }}\\&=\ SEZ {g _ {I, j + 1} U _ {I, j + 1} + G _ {IJ} U _ {I, j-1} + G _ {I + 1, J} U _ {I + 1, J} + G _ {IJ} U _ {I-1, j }}\\&-\ sex {g _ {I, j + 1} + G _ {IJ} + G _ {I + 1, j} + G _ {I, j} U _ {IJ}, \ EEx $. To ensure the computing stability, we use a semi-implicit solution to solve the problem: $ \ bee \ label {2: Q} \ Bea Q (U ^ {n + 1 }_{ IJ}) & =\ SEZ {G ^ N _ {I, J + 1} U ^ N _ {I, j + 1} + G ^ N _ {IJ} U ^ N _ {I, j-1} + G ^ N _ {I + 1, J} U ^ N _ {I + 1, J} + G ^ N _ {IJ} U ^ N _ {I-1, j }}\\& + \ sex {G ^ N _ {I, j + 1} + G ^ N _ {IJ} + G ^ N _ {I + 1, j} + G ^ N _ {I, j} U ^ {n + 1 }_{ IJ }, \ EEA \ EEE $ where $ G ^ N _ {IJ} $ is selected as follows (dependent on its corresponding $ U _ {IJ} $, one Direction is backward or forward, and the other direction is center difference ):

(1) $ U ^ N _ {I, j + 1} $ coefficient $ \ Bex C_1 \ equiv g _ {I, J + 1 }=\ frac {1} {\ SQRT {\ sex {u _ {I, J + 1}-U _ {IJ} ^ 2 + \ sex {\ frac {u _ {I + 1, J + 1}-U _ {I-1, J + 1 }}{ 2 }}^ 2 }}; \ EEx $

(2) $ U ^ N _ {I, J-1} $ coefficient $ \ Bex C_1 \ equiv g _ {I, J + 1 }=\ frac {1} {\ SQRT {\ sex {u _ {IJ}-U _ {I, j-1} ^ 2 + \ sex {\ frac {u _ {I + 1, J-1}-U _ {I-1} {2 }}; \ EEx $

(3) $ U ^ N _ {I + 1, J} $ coefficient $ \ Bex C_1 \ equiv g _ {I, J + 1 }=\ frac {1} {\ SQRT {\ sex {u _ {I + 1, j}-U _ {IJ} ^ 2 + \ sex {\ frac {u _ {I + 1, J + 1}-U _ {I + 1, j-1 }}{ 2 }}^ 2 }}; \ EEx $

(4) $ U ^ N _ {I-1, j} $ coefficient $ \ Bex C_1 \ equiv g _ {I, J + 1 }=\ frac {1 }{\ SQRT {\ sex {u _ {IJ}-U _ {I-1, j }}^ 2 + \ sex {\ frac {u _ {I-1, J + 1}-U _ {I-1} {2 }}. \ EEx $ substitute \ eqref {2: Q} into \ eqref {2: Scheme }, $ \ Bex & \ SEZ {1 + \ MU \ Tau \ Delta _ \ ve (U ^ N _ {IJ })} U ^ {n + 1 }_{ IJ} = u ^ N _ {IJ} + \ Tau \ Delta _ \ ve (U ^ N _ {IJ }) \ cdot \ & \ quad \ cdot \ SEZ {\ MU \ sex {c_1u ^ N _ {I, j + 1} + c_2u ^ N _ {I, j-1} + c_3u ^ N _ {I + 1, J} + c_4u ^ N _ {I-1, j} + (I _ {IJ}-C ^ N_1) ^ 2-(I _ {IJ}-C ^ N_2) ^ 2}, \ EEx $ where $ C ^ N_1 $, $ C ^ N_2 $ by \ eqref {1: c1c2} is discretization: $ \ Bex C ^ N_1 = \ frac {\ sum _ {IJ} I _ {IJ} \ SEZ {1-H _ \ ve \ sex {u ^ N _ {IJ }} }}{ \ sum _ {IJ} \ SEZ {1-H _ \ ve \ sex {u ^ N _ {IJ }}}}, \ quad C ^ N_2 = \ frac {\ sum _ {IJ} I _ {IJ} H _ \ ve \ sex {u ^ N _ {IJ }}{\ sum _ {IJ} H _ \ ve \ sex {u ^ N _ {IJ }}}. \ EEx $

 

3. Numerical Test

Select $ \ ve = 1.0 $, $ \ mu = 250 $, $ \ Tau = 0.1 $, and use the following Matlab code: Clear all;

Close all;

CLC;

 

Img‑imread('brain.bmp ');

IMG = double (rgb2gray (IMG ));

 

Figure (1 );

Imshow (uint8 (IMG ));

 

[NY, nx] = size (IMG); % get image size

 

% Set the initial curve to a circle

% Center of the initial circle

C_ I = floor (NY/2 );

C_j = floor (NX/2 );

% Radius of the initial circle

R = C_ I/3;

 

% Initialize U as a distance function

U = zeros ([NY, nx]);

For I = 1: NY

For j = 1: NX

U (I, j) = r-SQRT (i-c_ I). ^ 2 + (j-c_j). ^ 2 );

End

End

 

% Place the initial circular curve on the original image

Figure (2 );

Imshow (uint8 (IMG ));

Hold on;

[C, H] = contour (u, [0 0], 'R ');

 

% Initialization parameters

Epsilon = 1.0; % Heaviside function parameter settings

MU = 250; % arc length weight

Dt = 0.1; % time step

Nn = 0; % initialize the number of output images

 

% Iteration start

For N = 1:1000

% Calculate the normalized Heaviside Function

H_u = 0.5 + 1/PI. * atan (U/epsilon );

% Calculated by the current U parameter C1, C2

C1 = sum (1-h_u). * IMG)/sum (1-h_u ));

C2 = sum (h_u. * IMG)/sum (h_u ));

% Update u with current C1, C2

Delta_epsilon = 1/pI * Epsilon/(PI ^ 2 + Epsilon ^ 2); % regularization of Delta functions

M = DT * delta_epsilon; % product of the temporary Parameter M storage time step and Delta Function

% Calculate the weight of the four neighboring points

C1 = 1./SQRT (EPS + (U (:, [2: NX, nx])-U). ^ 2

+ 0.25 * (U ([2: NY, NY], :)-U ([: ny-1], :)). ^ 2 );

C2 = 1./SQRT (EPS + (u-u (:, [: nx-1]). ^ 2

+ 0.25 * (U ([2: NY, NY], [: nx-1])-U ([: ny-1], [: nx-1]). ^ 2 );

C3 = 1./SQRT (EPS + (U ([2: NY, NY], :)-U). ^ 2

+ 0.25 * (U ([2: NY, NY], [2: NX, nx])-U ([2: NY, NY], [1, 1: nx-1]). ^ 2 );

C4 = 1./SQRT (EPS + (u-u ([: ny-1], :)). ^ 2

+ 0.25 * (U ([: ny-1], [2: NX, nx])-U ([: ny-1], [: nx-1]). ^ 2 );

C = 1 + Mu * m. * (C1 + C2 + C3 + C4 );

U = (U + M *

(MU * (C1. * U (:, [2: NX, nx]) + C2. * U (:, [: nx-1])

+ C3. * U ([2: NY, NY], :) + C4. * U ([: ny-1], :))

+ (Img-c1). ^ 2-(Img-c2). ^ 2)

)./C;

% The curve and shard constant curve are displayed every two hundred running times.

If Mod (n, 200) = 0

Nn = nn + 1;

F = IMG;

F (u> 0) = C1;

F (u <0) = c2;

Figure (Nn + 2 );

Subplot (1, 2); imshow (uint8 (f ));

Subplot (1, 2); imshow (uint8 (IMG ));

Hold on;

[C, H] = contour (u, [0, 0], 'R ');

Hold off;

End

End

To achieve the division of human brain slices, see:

 

Note:

1. For ease of writing, some parts of the MATLAB program have been separated, such as $ C_1, \ cdots, C_4 $, and $ U $. You cannot segment it during writing.

2. From the above procedures, we can see that multi-purpose matrix operations in Matlab are much faster than pure iterative operations.

 

4. Thank you

The author wowould like to thank Dr. s.j. li and dr. p. li for suggesting the book by D. k. wang and etc. and special thank goes to dr. p. li for exploiting this MATLAB procedure.