6. Code Matlab:

exemple permettant d'éxécuter dans l'environnement Matlab (6.5) le programme matlab de l'application :

1- Lecture de l'image en spécifiant le chemin complet vers son fichier,
par exemple: c:/Bechar_Project/gvf4_fichiers/brain2.gif

 >> w=get_gray_level_image_matrix('c:/Bechar_Project/gvf4_fichiers/brain2.gif');

2- Spécification des paramètres et initilisation du champ GVF: Par exemple:
 >> mu=0.25; taw=20; iterMax=20; epsil=0.0001;
 >> N=size(w);
 >> u_0=rand(256,256); v_0=rand(256,256);

3- Calcul et affichage du champ GVF :

 >> [u,v,V,iterations,precision]=gvf_field(mu,taw,u_0,v_0,w,iterMax,epsil);
 
4- Spécifier avec la souris (par un seul clic à la fois sur l'image) des points
du snake initial:

 >> [x0,y0]=specify_initial_snake_points(w);

5- Spécification des paramètres de GVF snake: Par exemple:

 >>  alpha=0.08;   beta=0;  gamma=50; taw=0.0001; iterMax=10; epsil=0.00001;

6- spécification des points fixés dans le snake: Par exemple:

 >> Fixed=[1,5,15,25,30];

(cela suppose bien sûr que vous avez sélectionné au moins 30 points dans le snake
initial).

5- Calcul et affichage du snake GVF final:

>> [x,y,iter,prec]=gvf_snake(alpha,beta,gamma,taw,x0,y0,Fixed,u,v,w,iterMax,epsil);

.

function w=get_gray_level_image_matrix(c)

% w=read_gray_level_image_matrix(c)
% returns the gray-level matrix of the image stocked in file c.
% c: path to the file where the image is stocked

[im,clmp]=imread(c);
w=ind2gray(im,clmp);

function [u,v]=iter_gvf_field(mu,taw,u_t,v_t,I)

mu_taw=mu*taw;
[N,M]=size(I);

gradIx=diff(I);
gradIx(N,:)=(2*gradIx(N-1,:)+gradIx(N-2,:))/3;

gradIy=diff(I');
gradIy(M,:)=(2*gradIy(M-1,:)+gradIy(M-2,:))/3;
gradIy=gradIy';

% second
% computation of f as the square of modulus of gradient of image I
f=(gradIx.^2+gradIy.^2);
 

%computation of the matrix fx (partial derivative of f with respect to x
fx=diff(f);
fx(N,:)=(2*fx(N-1,:)+fx(N-2,:))/3;

%computation of the matrix fy (partial derivative of f with respect to y
fy=diff(f');
fy(M,:)=(2*fy(M-1,:)+fy(M-2,:))/3;
fy=fy';

%computation of the second member in the equation system of the GVF field
Fu=2*taw*(fx.^2+fy.^2).*(u_t-fx)-2*u_t;
Fv=2*taw*(fx.^2+fy.^2).*(v_t-fy)-2*v_t;

u=inversion(N,M,mu_taw,Fu);
u=real(u); % suppression of the imaginay part of u due to calculus approximations
v=inversion(N,M,mu_taw,Fv);
v=real(v); % suppression of the imaginay part of u due to calculus approximations
return

% computation of u at iteration t

function u=inversion(N,M,mu_taw,F)
alpha=4*mu_taw;
beta=-(2+8*mu_taw);

% computation of the fast Fourier transorm of F
fftF=fft2(F);

% computation of the fast Fourier transorm of u
for k=1:N,
for j=1:M,
a(k,j)=alpha*(cos(2*pi*(k-1)/N)+cos(2*pi*(j-1)/M))+beta;
end
end

fftu=fftF./a;

u=ifft2(fftu);

function [u,v,V,iterations,precision]=gvf_field(mu,taw,u_0,v_0,I,iter,epsil)

% [u,v,V,iter,prec]=gvf_field(mu,taw,u_t,v_t,I,iter,epsil): returns the GVF field
% mu: mu parameter of GVF Algorithm
% iter: maximum number of iterations before reaching the steady point of
% GVF field
% epsil: precision: indicates that the steady point of GVF field has been reached
% u_0: initial value of u ; the first field component
% v_0: initial value of v; the second field component
% I : gray-level image
% sigma: gaussian deviation parameter used to compute the edge map using
% the formula : f=-Eext, where Eext is the gradient of convlolution of image I with the gaussian
% of deviation sigma

iterations=0;
precision=100000;
% computation of the edge map
u=0;
v=0;
while((iterations<iter)&(precision>epsil)),
[u,v]=iter_gvf_field(mu,taw,u_0,v_0,I);
precision=norm([u-u_0,v-v_0]);
u_0=u;
v_0=v;
iterations=iterations+1,

end
V=sqrt(u.^2+v.^2);

figure;
imshow(V);

return

function [x0,y0]=specify_initial_snake_points(w)

%[x0,y0]=specify_initial_snake_points(w)
% allows for the user to specify the points to be part of the initial snake
% in an image given as a matrix of gray levels
% it retruns their coordiantes x0 and y0
[M,N]=size(w);
i=1;
while (1)
[c,b]=imcrop(w);
if (b(1)<0)|(b(1)>M)|(b(2)<0)|(b(2)>N)
break;
else
y0(i)=b(1);
x0(i)=b(2);
i=i+1;
end

end

function [X,Y]=curv_param(x0,y0)

x00=x0;
y00=y0;

N=length(x0);
L(1)=0;

X(1)=x0(1);
Y(1)=y0(1);
Ns=1; % the number of already computed points is one

for i=2:N,
if x00(i)==x00(i-1),
L(i)=L(i-1)+abs(y00(i)-y00(i-1));
ns=integer(L(i)-Ns);
for j=1:ns,
X(Ns+j)=x00(i);
Y(Ns+j)=y00(i)+j;
end
 

else
L(i)=L(i-1)+sqrt((x00(i)-x00(i-1))^2+(y00(i)-y00(i-1))^2);
ns=integer(L(i)-Ns); % the number of points to be computed at iteration i
v=1+(y00(i)-y00(i-1))^2/(x00(i)-x00(i-1))^2;
v=sqrt(v);
if x00(i)>x00(i-1),
for j=1:ns,
X(Ns+j)=x00(i-1)+(Ns+j-L(i-1))/v;
Y(Ns+j)=(y00(i)-y00(i-1))*(X(Ns+j)-x00(i-1))/(x00(i)-x00(i-1))+y00(i-1);
end
else
for j=1:ns,
X(Ns+j)=x00(i-1)-(Ns+j-L(i-1))/v;
Y(Ns+j)=(y00(i)-y00(i-1))*(X(Ns+j)-x00(i-1))/(x00(i)-x00(i-1))+y00(i-1);
end
end
end

Ns=Ns+ns; % Ns: the number of already computed points until iteration i, also equal to the
% integer value of the curve length at iteration i

end

%figure;
%plot(X,Y);

function [x,y]=iter_gvf_snake(alpha,beta,gamma,taw,x0,y0,Xfixed,Yfixed,u,v)

L=length(x0);
for h=1:L,
q=x0(h);
r=y0(h);
u_i(h)=-gamma*taw*interpol(u,q,r)-q;
v_i(h)=-gamma*taw*interpol(v,q,r)-r;
end

%computation of the fft of the second member
tfui=fft(u_i);
tfvi=fft(v_i);

for m=1:L,
a(m)=2*taw*beta*cos(4*pi*(m-1)/L)+2*taw*(-4*beta+alpha)*cos(2*(m-1)*pi/L)+taw*(6*beta-2*alpha)-1;
end

% - computation of the fft of the sname
tfx=tfui./a;
tfy=tfvi./a;

% computation of the sname using the ifft
x=real(ifft(tfx));
y=real(ifft(tfy));
t=length(Xfixed);

for j=1:L,
replaced(j)=0;
end

for i=1:t,
min=0; dmin=100000000000;
for j=1:L,
if ~replaced(j),
d=distance(Xfixed(i),Yfixed(i),x(j),y(j));
if d<dmin,
min=j; dmin=d;
end
end

end
x(min)=Xfixed(i);
y(min)=Yfixed(i);
end

return
 

%interpolation of gvf field :(u,v)
function z=interpol(u,x,y)
[rangeX,rangeY]=size(u);
if (x>rangeX)|(x<1)|(y>rangeY)|(y<1),
z=0;
return
end

x1=integer(x);
y1=integer(y);
if (x1<=1)|(x1>=rangeX),
dux1=0;
else
dux1=u(x1+1,y1)-u(x1,y1);
end
if (y1<=1)|(y1>=rangeY),
duy1=0;
else
duy1=u(x1,y1+1)-u(x1,y1);
end
z=u(x1,y1)+(x-x1)*dux1+(y-y1)*duy1;

function d=distance(x1,y1,x2,y2)
d=sqrt((x1-x2)^2+(y1-y2)^2);

function [x,y,iterations,prec]=gvf_snake(alpha,beta,gamma,taw,x0,y0,Fixed,u,v,I,iterMax,epsil)

%[x,y,iter,epsil]=gvf_snake(alpha,beta,taw,x_0,y_0,u,v,I,iterMax,epsil)
% alpha : alpha parameter of the GVF snake model
% beta : beta parameter of the GVF snake model
% gamma: scale parameter
% taw : time pace parameter
% x0 and y0 : initial snake
% Fixed : vector of the indices in the initial snake pointys of the fixed
% points
% V=[u,v]: GVF field precedently calculated
% I : matrix of a gray-level image
% iterMax: superior bound of the number of iterations
% epsil: wanted precision

fxd=length(Fixed);

for i=1:fxd,
Xfixed(i)=x0(Fixed(i));
Yfixed(i)=y0(Fixed(i));
end

[rangeX,rangeY]=size(u);
prec=1000000;
iterations=0;
xt=x0;
yt=y0;
xt1=x0;
yt1=y0;

[xt,yt]=curv_param(xt1,yt1);
L=length(xt),
while((iterations<iterMax)&(prec>epsil)),
[xt,yt]=curv_param(xt1,yt1);
[xt1,yt1]=iter_gvf_snake(alpha,beta,gamma,taw,xt,yt,Xfixed,Yfixed,u,v);
prec=norm([xt1-xt,yt1-yt]);
xt=xt1; yt=yt1;
L=length(xt);
iterations=iterations+1,
end

for i=1:L,

x1(i)=integer(xt(i));
if xt(i)-x1(i)>=0.5,
x1(i)=x1(i)+1;
end
xt2(i)=max(x1(i),1);
xt2(i)=min(xt2(i),rangeX);

y1(i)=integer(yt(i));
if yt(i)-y1(i)>=0.5,
y1(i)=y1(i)+1;
end
yt2(i)=max(y1(i),1);
yt2(i)=min(yt2(i),rangeY);
end

g=I; f=sqrt(u.^2+v.^2);
for i=1:L,
g(xt2(i),yt2(i))=1;
f(xt2(i),yt2(i))=1;
end

figure;
imshow(f);
figure;
imshow(g);
x=xt2;
y=yt2;

return

  Page précédente

  Première Page