%load sdf
load shape

dt = 0.2;
max_it = 10000;
eps=1e-12;

phi = shape;
for i = 1:max_it
   
    %compute first and second central differences
    phi_cx = 0.5*(circshift(phi, [-1 0]) - circshift(phi, [1 0]));
    %also phi_cx = imfilter(phi, [-0.5; 0.0; 0.5], 'circular');
    
    phi_cy = 0.5*(circshift(phi, [0 -1])-circshift(phi, [0 1]));
    phi_cxx = circshift(phi, [1 0])- 2.*phi + circshift(phi, [-1 0]);
    phi_cyy = circshift(phi, [0 1])- 2.*phi + circshift(phi, [0 -1]);
    phi_cxy = 0.5*(circshift(phi_cx, [0 -1])-circshift(phi_cx, [0 1]));
    
    gradMag = (phi_cx.^2 + phi_cy.^2).^0.5;
    
    numer = -(phi_cxx.*phi_cy.^2 - 2*phi_cx.*phi_cy.*phi_cxy + phi_cyy.*phi_cx.^2);
    denom = ((phi_cx.^2 + phi_cy.^2).^1.5 + eps);
    curvature = numer./denom;
    
    phi = phi - dt*curvature.*gradMag;
         
    if mod(i,50)==0
        % signed distance reinitialization
        phi = double((phi > 0).*(bwdist(phi < 0)-0.5) - (phi < 0).*(bwdist(phi > 0)-0.5));
    end
    
    if mod(i,100)==0
        imagesc(phi>0);
        pause(0.01);
    end
        
end