$ \ 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.