% To run this script, just type levelset at the >> prompt. 
x=double(imread('simpleimage.tif'));
figure(1)
imagesc(x);
axis('image');
colormap('gray');
maxx=max(max(x));
minx=min(min(x));
x=(x-minx)./(maxx-minx);
[m,n]=size(x);
sign=0.01;
xn=x+sign*randn(m,n);
filsiz=9;
flup=ceil(filsiz/2);
sigma=0.8;
for i=1:filsiz
    for j=1:filsiz
        gauss(i,j)=exp(-((i-flup)^2+(j-flup)^2)/(2*sigma^2));
    end
end
figure(2)
imagesc(xn);
axis('image');
colormap('gray');
y=conv2(xn,gauss);
figure(3)
imagesc(y);
axis('image');
colormap('gray');
% Gradient calculation
[m,n]=size(y);
gyx=zeros(m,n);
for i=1:n
    for j=1:n-1
        gyx(i,j)=y(i,j+1)-y(i,j);
    end
end
figure(4)
imagesc(gyx);
axis('image');
colormap('gray');
gyy=zeros(m,n);
for i=1:n-1
    for j=1:n
        gyy(i,j)=y(i+1,j)-y(i,j);
    end
end
figure(5)
imagesc(gyy)
axis('image');
colormap('gray');
delta=0.01;
speedfunc=1./(1+(gyx.^2+gyy.^2)/delta);
figure(6)
imagesc(speedfunc);
axis('image');
colormap('gray');
% Initialization of S
for i=1:m
    for j=1:n
        D=55;
        if ((i-m/2)^2 + (j-n/2)^2) < D^2
            S(i,j)=(-D^2+(i-m/2)^2+(j-n/2)^2)/D^2;
            %S(i,j)=-1;
        else
            S(i,j)=(((i-m/2)^2+(j-n/2)^2)-D^2)/D^2;
            %S(i,j)=1;
        end
    end
end
Sini=S;
figure(7)
imagesc(S)
axis('image');
colormap('hsv');
% Start iterating
Sx=zeros(m,n);
Sy=zeros(m,n);
for k=1:7
    k
    for i=1:m
        for j=2:n-1
            Sx(i,j)=S(i,j+1)-S(i,j-1);
        end
    end
    for i=2:m-1
        for j=1:n
            Sy(i,j)=S(i+1,j)-S(i-1,j);
        end
    end
    alpha=1;
    for i=1:m
        for j=1:n
            gradS(i,j)=Sx(i,j)^2+Sy(i,j)^2;
            S(i,j)=S(i,j)-alpha*speedfunc(i,j)*sqrt(Sx(i,j)^2+Sy(i,j)^2);
        end
    end
    S=S./max(max(S));
    figure(9)
    imagesc(gradS);
    axis('image');
    colormap('gray');
    figure(10)
    quiver(Sx,Sy);
    axis('image');
    figure(8)
    imagesc(S);
    axis('image');
    colormap('hsv');
    pause(0.3)
end

    



