function test_ac
im=imread('flower.tif');
figure
imshow(im);
x=130+75*cos(0:0.1:2*pi)';%initial snake x,y are column vectors
y=130+75*sin(0:0.1:2*pi)';
hold on 
plot([x;x(1)],[y;y(1)],'g')
hold on

Nx=length(x);
Ny=length(y);

mu=1;
alpha=0.001;
beta=0.37;
dt=2;
iterations=5000;

a=dt*(2*alpha+6*beta)+1;
b=dt*(-alpha-4*beta);
c=dt*beta;
P=diag(repmat(a,1,Nx));%P=id-dt*A
P=P+diag(repmat(b,1,Nx-1),1)+diag(b,-Nx+1);
P=P+diag(repmat(b,1,Nx-1),-1)+diag(b,Nx-1);
P=P+diag(repmat(c,1,Nx-2),2)+diag([c,c],-Nx+2);
P=P+diag(repmat(c,1,Nx-2),-2)+diag([c,c],Nx-2);
P=inv(P);

H=fspecial('gaussian',[10,10],10);
im=imfilter(im,H);
[FX,FY]=gradient(double(im));%FX,FY are Nx x Ny matrices
%H=fspecial('gaussian',[10,10],30);
%FX=imfilter(FX,H);
%FY=imfilter(FY,H);
[FX,FY]=gradient((FX.^2+FY.^2).^(1/2));

Ex=zeros(Nx,1);
Ey=zeros(Nx,1);
for i=1:iterations
    FX_interpol=griddedInterpolant(FX,'linear');
    FY_interpol=griddedInterpolant(FY,'linear');
    for j=1:Nx
        Ex(j)=FX_interpol(y(j),x(j));
        Ey(j)=FY_interpol(y(j),x(j));
    end
    x=P*(x+dt*mu*Ex);
    y=P*(y+dt*mu*Ey);
    if mod(i,5)==0
        plot([x;x(1)],[y;y(1)],'b')
    end
end
plot([x;x(1)],[y;y(1)],'r')
