M = csvread('hilbert.csv');
x = abs(M(:,1)');
y = abs(M(:,2)');
% calculate intermediate points so that we get an approximate equal spaced
% line of points
n = length(x);
xi = []; yi = [];

stepsize = 0.3;%new stepsize

original_stepsize = sqrt((x-[x(2:end),x(1)]).^2+(y-[y(2:end),y(1)]).^2);%os(i) = d(i,i+1)
for i=1:n
    t = linspace(0,1,original_stepsize(i)/stepsize);
    t = t(1:end-1);
    if i ~= n
        xi = [xi, x(i)*(1-t)+x(i+1)*t];
        yi = [yi, y(i)*(1-t)+y(i+1)*t];
    else %i==n
        xi = [xi, x(n)*(1-t)+x(1)*t];
        yi = [yi, y(n)*(1-t)+y(1)*t];
    end
end
%% animation
h = figure(1);
set(h, 'Position', [100, 100, 600, 600]);
filename = 'curve_shortening.gif';
flag_first = 1;

timestepsize = 0.5;
totalit = 0;%total nubmer of iterations
for k=1:300

plot(xi,yi,'-')
text(50,275,[num2str(totalit),' ',num2str(k)])
axis([-10,265,-10,265]);

%calculate first and second derivatives

    for l=1:k^0.8
        ddxi = [xi(2:end),xi(1)]-2*xi+[xi(end),xi(1:end-1)];
        ddyi = [yi(2:end),yi(1)]-2*yi+[yi(end),yi(1:end-1)];
        xi = xi+timestepsize*ddxi;
        yi = yi+timestepsize*ddyi;
        totalit = totalit+1;

        if mod(k,1) == 0
            r = randi(length(xi),1,round(length(xi)*0.0005));%remove 0.05percent of points
            xi(r) = [];
            yi(r) = [];
        end

    end
%capture image
    if mod(k,2) == 1
        set(gcf,'color',[1 1 1]);%white background
        axis off
        drawnow
        frame = getframe(1);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);
        if flag_first == 1;
            imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',0.05);
            flag_first = 0;
        else
            imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.05);
        end
    end

%pause(0.1)

end