% clear all, clc;

% parameter setting
sigma = 0;   % gaussianBlur
GradientOn = 1;  % 1 : Gradient on  0 : Gradient off

GVFOn = 1;   % 1 : GVF  0 : SVF
Num = 150;      % number of GVF iterations
mu = 0.1;    % parameter for GVF (more noise, increase mu)

iter = 50;   % number of Snake iterations 
alpha = 0.05;   % internal weight that control snake's tension
beta = 0;     % internal weight that control snake's rigidity
gamma = 1;    % step size in one iteration
kappa = 0.5;   % external force weight
Dmin = 0;   % min raster of the snake
Dmax = 1;   % max raster of the snake

Ra_min = 20;
Ra = 1000;   % maximum area
Rc = 0.9;   % convexity threshlod
Rl = 2;    % threshold ratio between length and width

se = strel('disk', 3);  % morphology structure element

r = 20;  % radium of initial circle contour
x0 = 80;
y0 = 40;   % centroid of initial contour

timer = 1;


d = 10;
order = 2;

II = imread('test8.jpg');
I = rgb2gray(II);
bw = im2bw(I, graythresh(I));
[s1, s2] = size(bw);
s_1 = [0, s1, s1, 0];
s_2 = [0, 0, s2, s2];  % vertex of image boundary

% blured and draw the picture
if sigma ~= 0 
   f = gaussianBlur(I,sigma);
else
   f = I;
end;

% graient on the image in order to get edges
if GradientOn 
	f2 = abs(gradient2(double(f)));
else 
    f2 = f;
end;

% vector field
if GVFOn                          
    % GVF field
    [u, v] = GVF(f2, mu, Num);   % gradient vector force
else
    % standard vector field
    [u, v] = gradient2(f2);   % calculates standard external force vector
                              % field using gradient
end

% normalize vectors in vector field
mag = sqrt(u.*u + v.*v);
px = u ./ (mag + 1e-10);
py = v ./ (mag + 1e-10);

xSpace = (1:size(bw,1)/64:size(bw,1));
ySpace = (1:size(bw,2)/64:size(bw,2));
qx = interp2(px, xSpace, ySpace');
qy = interp2(py, xSpace, ySpace');
quiver(xSpace, ySpace, qx, qy); 
axis('ij'); axis equal; axis off;

% snake deformation
% imshow(I);hold on,
t = 0:0.05:6.28;
x = x0 + r*cos(t);
y = y0 + r*sin(t);
[x, y] = polybool('intersection', s_2, s_1, x, y);
                             % in case polygon vertex is outside of image
                             % x: horizontal, y: vertical
% plot(x,y,'r');

figure,imshow(II);
x = x';
y = y';
[x,y] = snakeinterp(x, y, Dmax, Dmin);
snakedisp(x,y,'r') 
pause(1);

for i=1 : ceil(iter / 5),
   if i <= floor(iter / 5) 
      [x, y] = snakedeform(x, y, alpha, beta, gamma, kappa, px, py, 5);
   else 
      [x, y] = snakedeform(x, y, alpha, beta, gamma, kappa, px, py, iter - floor(iter / 5)*5);
   end;
   [x, y] = snakeinterp(x, y, Dmax, Dmin);
   snakedisp(x,y,'y')
   pause(0.1);
end
hold on,
plot(x,y,'g');

% bw = roipoly(I,x,y);
% xx = ceil(x); yy = ceil(y);
% f3 = I(xx(:), yy(:)) == 0;
% imshow(I),hold on,plot(xx,yy,'r.')
bw = im2bw(I, graythresh(I));
[v, h] = size(bw);
xx = ceil(x); yy = ceil(y);
len = length(xx);
for i = 1:len
    if xx(i) <= h && yy(i) <= v
        bw(yy(i),xx(i)) = 0;
    end
end
% figure,imshow(bw);


