Contour actif

clear;

% I=double(rgb2gray(imread('Synth.jpg')));
% I=double(rgb2gray(imread('Synth2.jpg')));
I=double(rgb2gray(imread('fruitfly.jpg')));
% I=double(rgb2gray(imread('objects.jpg')));

imshow(uint8(I));
rect=getrect(gcf);
taille=(size(I));
xmin=rect(1);
ymin=rect(2);
width=rect(3);
height=rect(4);
% [xmin ymin width height]=GETRECT((I))

NBpts=120;

%% placement des points du contours original
NBpts=floor(NBpts/4);
ecart_largeur=floor((width)/NBpts);
ecart_hauteur=floor((height)/NBpts);
NBpts=NBpts*4;
for (i=1:NBpts)
if (i<=NBpts/4) %% horizontale haut
p(i,2)=round((i-1)*ecart_largeur+xmin);
p(i,1)=round(ymin);
end
if ((i>NBpts/4)&(i<=2*NBpts/4)) %% verticale droite
p(i,2)=round(xmin+width);
p(i,1)=round(((i-1)-(NBpts/4))*ecart_hauteur+ymin);
end
if ((i>2*NBpts/4)&(i<=3*NBpts/4)) %% Horizontale basse
p(i,2)=round(xmin+width-((i-1)-(2*NBpts/4))*ecart_largeur);
p(i,1)=round(ymin+height);
end
if (i>3*NBpts/4) %% Verticale gauche
p(i,2)=round(xmin);
p(i,1)=round(ymin+height-((i-1)-(3*NBpts/4))*ecart_hauteur);
end
end

%% affichage des pts
Itmp=I;
for i=1:NBpts
Itmp(p(i,1),p(i,2))=0;
end
imshow(uint8(Itmp))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Calcul algo glouton

%image gradiant
Igrad=-(filter2(fspecial('laplacian'),I));

Eelas=zeros(3,3);
Erigid=zeros(3,3);
for cpt=1:1000
for i=1:NBpts
M=9;
m=1;
for a=-1:1
for b=-1:1

%calcul energie elastique (interne)
if(i==1)
Eelas(a+2,b+2)=(p(i,1)+a-p(NBpts,1))^2+(p(i,2)+b-p(NBpts,2))^2;
else
Eelas(a+2,b+2)=(p(i,1)+a-p(i-1,1))^2+(p(i,2)+b-p(i-1,2))^2;
end

%% calcul energie rigidite (interne)
if(i==1)
Erigid(a+2,b+2)=(p(NBpts,1)-2*(p(i,1)+a)+p(i+1,1))^2+(p(NBpts,2)-2*(p(i,2)+b)+p(i+1,2))^2;
else if (i==NBpts)
Erigid(a+2,b+2)=(p(i-1,1)-2*(p(i,1)+a)+p(1,1))^2 + (p(i-1,2)-2*(p(i,2)+b)+p(1,2))^2;
else
Erigid(a+2,b+2)=(p(i-1,1)-2*(p(i,1)+a)+p(i+1,1))^2+(p(i-1,2)-2*(p(i,2)+b)+p(i+1,2))^2;
end
end
m=m+1;
end
end

%calcul energie externe (gradiant)
Egrad(1:3,1:3)=Igrad(p(i,1)-1:p(i,1)+1,p(i,2)-1:p(i,2)+1);

E=0.4*Erigid+1.4*Eelas+3*Egrad;

[x,y]=find(E==min(min(E)));% recherche du minimum de tout les energies pour donnée les nouvelles valeurs
% a p.

p(i,1)=x(1)+p(i,1)-2;%-2 car decalage de la matrice (les indices)
p(i,2)=y(1)+p(i,2)-2;%-2 car decalage de la matrice (les indices)
end
if((cpt/20)==round(cpt/20))
%affichage des pts
Itmp=I;
Itmp(:,:,2)=I;
Itmp(:,:,3)=I;
i=1;
Itmp(p(i,1),p(i,2),1)=0;
Itmp(p(i,1),p(i,2),2)=255;
Itmp(p(i,1),p(i,2),3)=0;
for i=2:NBpts
Itmp(p(i,1),p(i,2),1)=255;
Itmp(p(i,1),p(i,2),2)=0;
Itmp(p(i,1),p(i,2),3)=0;
end
imshow(uint8(Itmp))
pause
end
end

voila il y a tout le code a toi de faire le trie