% This snake program makes use of sparse linear solver
% a simple trick leading to some approximation is used to reduce the bandwidth
% of the coefficient matrix to only a tridiagonal one
close all;
clear all;
clc;

h=256;
w=256;
I = zeros(h,w);
I(round(h/2),round(w/2))=1;
I = double(bwdist(I)<round(h/4));
I = imfilter(I,fspecial('gaussian',[17 17],3),'same','conv','replicate');
figure(1),imagesc(I),colormap(gray),axis image;
set(gcf,'doublebuffer','on');

N=500; % number of snake points
alpha=1;
tstep=1;
N_iter=500;
f=50;

[Ix,Iy] = gradient(I);
S = -f*(Ix.*Ix + Iy.*Iy);
[Sx,Sy] = gradient(S);
Smag = sqrt(Sx.^2 + Sy.^2)+eps;
Sx(:) = Sx./Smag;
Sy(:) = Sy./Smag;

theta = linspace(0,2*pi,N);
theta = theta(:);
x = w/2 + 10 + (h/3)*cos(theta);
y = h/2 - 10 + (h/4)*sin(theta);

figure(1), hold on;
plot(x,y,'r');

D=[-tstep*alpha*ones(N,1) (1+2*tstep*alpha)*ones(N,1) -tstep*alpha*ones(N,1)];
D(2,3)=D(2,3)-tstep*alpha;
D(end-1,1)=D(end-1,1)-tstep*alpha;

A = spdiags(D,[-1 0 1],N,N); % A is a tri diagonal matrix - O(n) algorithm to solve the system
bx = zeros(N,1);
by = zeros(N,1);
tic;
for iter=1:N_iter,
    bx(:) = x - tstep*interp2(Sx,x,y);
    by(:) = y - tstep*interp2(Sy,x,y);
    
    % solve by linear solver - inbuilt Matlab solver
    x = A\bx;
    y = A\by;
    
    % Thomas algorithm will solve this in O(n)
%     x = thomas(A,bx);
%     y = thomas(A,by);
    
    % a circular shift to average out the effect of approximation
    % used in reducing the bandwidth of the coefficient matrix
    sft = round(N*rand);
    x(:) = circshift(x,sft);
    y(:) = circshift(y,sft);
    if rem(iter,10)==0,
        plot(x,y,'r'); drawnow;
    end
end
toc;
hold off;
figure(2),imagesc(I),colormap(gray),axis image;
hold on; plot(x,y,'g','Linewidth',2); hold off; drawnow;