snake1.m

Active Contours model(snake)

function snake1()
clc
clear all;close all;
ACCx=[];
ACCy=[];
alpha = .5; %string forces
beta = .5; %rod forces
gamma =75; %external forces
delta =100; %balloon forces(used to 00inflate the snake instead of contracting it,usefull in the pigheart image)
tao =0.02; %stepsize(step in which snake moves
w=1; 
InputPoints=1;
Nit =10; %# of iterations
I= imread('1.jpg');
[rows cols dims] = size(I);
if dims==3
    I = rgb2gray(I);
end

figure(1),imshow(I);
g=ginput(InputPoints);%argument defines the number of input points ,inja g yek matris ke tedad %satrha neshon dahndeh tedad noghate vorodiye
numOfPoints=size(g,1);% bode yek ya tedad noghati ke taiin shodeh ast ra moshakhas mikonad

Gauss=fspecial('gaussian', 5, 5); 
im = filter2(Gauss, I, 'same'); %convolving the image with the gaussian filter to get ppal structures
%im = conv2(I, Gauss,'same');
%im= imfilter(I, Gauss,'replicate'); conv2(inpImage,G,'same');
figure(2), imagesc(I),colormap gray %smoothed image
title('after Gauss filter');
[gradientX, gradientY] = gradient(im); % gradients in each direction
gradientXY=[gradientX, gradientY] ; 
figure(3), imagesc(gradientXY),colormap gray
title('image gradiant in x and y axis');


Potential = -w*sqrt(gradientX.*gradientX + gradientY.*gradientY); % potential energy
figure(4), imshow(Potential,[]);
title('image potential');
[FpotX,FpotY] = gradient(Potential); % assigned to each direction
%figure(6);imshow([FpotX,FpotY]);title('potential before negative');

% since it is not possible assign (-)to the last expression,because it will produce a single output
FpotX=-FpotX; % the potential force in X 
FpotY=-FpotY; % the potential force in Y 

g=sortrows(g);%ordering the points
numOfPoints=size(g,1);

%Method to initialize every snake. 
% initializing the perimeter of the snake S
InitSnakeradius=8;
NumberPointsPerimeter=20;
for k=1:numOfPoints;
    c=g(k,:);%center of circle
    % arguments(circle centre,number of points of snake,radius of circle)
    S=CircleSnake(c,NumberPointsPerimeter,InitSnakeradius); % S=perimeter(initial curve)

%Creating X an Y to mantain the configurations of the snake for all its
%states,initialized with S.
[u v]=size(S);
X(:,1)=S(:,1);
Y(:,1)=S(:,2);

%Iterations to compute the changes in the snake
fprintf('Starting iterations.......\n');
    for i=1:Nit;
        %b) External forces
        xy=[X(:,i),Y(:,i)]; % The snake corresponding to every iteration
        [p q]= size(xy);

        %second and fourth derivative by finite differences
        %since we will realize the same operation in all the iterations,we put
        %every position to use in the derivatives, in separate matrices.Otherwise
        %we will have to generate additional loops to control de positions

        xyplus1=[xy(2:p,:);xy(1,:)];
        xyplus2=[xy(3:p,:);xy(1:2,:)];
        xyminus1=[xy(p,:);xy(1:p-1,:)];
        xyminus2=[xy(p-1:p,:);xy(1:p-2,:)];

        %Computing internal forces
        % string forces
        Fstring=xyminus1-2*xy+xyplus1; 

        % rod forces
        Frod=-xyminus2+4*xyminus1-6*xy+4*xyplus1-xyplus2; 

        %Computing external forces(using the potential energy)
        %finding the corresponding potencial force,for the points of the snake(xy)
        %into the matrices FpotX,FpotY
        roundxy= round(xy); %because the index for FpotX or FpotY must be integers
            for k=1:p
                M(k)=FpotX(roundxy(k,2),roundxy(k,1));
                N(k)=FpotY(roundxy(k,2),roundxy(k,1));
            end;
        Fext=[M' N'];%natijeh ye matris misheh ke bordar aval(M) sotone aval on va bordar dovom sotone dovam on mishe
        % maghdar rang har kodam az 8 noghteh snake ra neghdari mikonad
        %balloon forces
         Fball=0; 
            if (delta ~=0)
                diff=(xyminus1-xyplus1);
                diffx=-diff(:,2);
                diffy=diff(:,1);
                numer=[diffx,diffy];
                denom=(((sum((numer.^2)'))').^0.5); % to avoid division by zero
                Fball=numer./[denom,denom];
            end;

        %the next position the snake is computed like the last position + the delta of moved distance, in x and y 
        NextXY=xy+tao*(alpha*Fstring+beta*Frod+gamma*Fext+delta*Fball); %formula Evolution of snake
        X(:,i+1)=NextXY(:,1); % X save all the Xvalues of the snake after each iteration
        Y(:,i+1)=NextXY(:,2);  % Y save all the Yvalues of the snake after each iteration
    end; % End of iterations

%Draw all snakes - emphasize first and last snake
if (k==1)
figure; imshow(I),axis image; hold on;
end

%showing the sequence of changing of the snake
for i=2:10:Nit; %i=1 will be the initial snake
    xy=[X(:,i),Y(:,i);X(1,i),Y(1,i)]; % the snake resulting of each iter
    line(xy(:,1), xy(:,2), 'Color','blue') ,hold on;
    drawnow;%draw the snake each iteration on the same image       
end;
% xy=[X(:,Nit+1),Y(:,Nit+1);X(1,Nit+1),Y(1,Nit+1)]; %final snake is Nit+1, since we contain in this matrix the initialSnake
% line(xy(:,1), xy(:,2), 'Color','yellow');hold on;
ACCx=[ACCx X];
ACCy=[ACCy Y];

end;
%Showing all the evolved snakes
figure(5),imshow(I),hold on;
plot(ACCx,ACCy,'Color','blue');