function [  ] = Homework_2_template_2011( )

    load abct.mat;
    I = abct;

    I = I - min(I(:));
    I = I./max(I(:));

    p = 20;     % number of points along curve
    r = 15;     % initial circle radius
    cen = [93,162];    % initial circle center

    s = [1:p]'; 
    x = cen(1) + r*cos(2*pi*s/p);
    y = cen(2) + r*sin(2*pi*s/p);

    imagesc(I);
    axis xy;
    colormap gray;
    hold on;
    plot([x', x(1)], [y', y(1)], 'LineWidth', 2);
    hold off

    for s=1:3 % subdivide 3 times during iteration for CS778
        [x, y] = EvolveSnake(I, x, y);
        [x, y] = Subdivide(x, y);
    end
end

function [x,y] = Subdivide(x, y)
%for CS 778 : fix this function. Perform midpoint subdivision on the input
%curve (x,y)
    S = speye(n,n);     % this S is w/o mid-point subdivision
    
    % below implementation is for S w/ mid-point subdivision
    S = zeros(2n,n);
    tmpDiagonal = zeros(n,1);tmpDiagonal(1)=1;
    for i = 1:(n-1)
        S = spdiags([],[],2n,n);    
    end
    
    x = S*x;
    y = S*y;
end

function [ x,y ] = EvolveSnake( image, x, y )
% runs 500 iterations of the snake evolution equation

    %You should be able to get good results with these values of a, b, dt,
    %but adjust if necessary.
    a = 0.02;   %alpha
    b = 0.005;  %beta
    dt = 0.2;
    gamma = 1.0/dt;

    h = fspecial('gaussian', 5, 1.5);
    Ismooth = imfilter(image, h);
    [gx, gy] = gradient(Ismooth);
    
    %pick a value for the exponent
    exponent = 2;
    
    %pick a value for the offset
    offset = 10;
    
    %pick a value for parameter K
    K = 2;
    
    % this will be the speed of the inflation
    f = 1./(1+ (sqrt(gx.^2 + gy.^2)./K).^exponent) - offset;
    
    n = numel(x);
    e = ones(n,1);
    
    x = reshape(x, [n,1]);
    y = reshape(y, [n,1]);
    
    I = speye(n,n);
    
    % for CS 578/778 : compute Df, Db, Dc, D2c, and A
    % be sure to impose periodic boundary conditions on Df, Db, Dc, and D2c
    
    Df = spdiags([-e, e, e],[0,1,-(n-1)], n, n);  % forward first difference matrix 
    Db = spdiags([e, -e, -e],[0,-1,(n-1)], n, n); % forward first difference matrix
    Dc = 0.5*spdiags([e, e, e],[0,1,-(n-1)], n, n);     % central first difference matrix
    D2c = spdiags([e, e, e],[0,-1,(n-1)], n, n);        % central second difference matrix
    % A should be a 5 by 5 sparse matrix with 5 non-zero diagonals
    %and it represents the smoothness of the curve
    A = spdiags([e,e,e,e,e],[0,1,2,-1,-2], n, n);
    
    [L, U] = lu(A+gamma*I);
    
    for i=1:500
        i
        
        dxds = Dc*x;
        dyds = Dc*y;
        len = sqrt(dxds.^2 + dyds.^2);
        
        %the unit tangent vector
        tx = dxds./len;
        ty = dyds./len;
        
        %the unit normal vector
        nx = -ty;
        ny = tx;
       
        %the external forces
        speed = ImageInterp(f, x, y);
        fx = speed.*nx; 
        fy = speed.*ny;
       
        tic;
        x = SolveLU(L,U,gamma*x - fx);
        y = SolveLU(L,U,gamma*y - fy);
        toc;
        
        imagesc(image);
        %imagesc(f);
        axis xy;
        colormap gray;
        hold on;
        plot([x', x(1)], [y', y(1)], '-+b', 'LineWidth', 2, 'MarkerEdgeColor', 'red');
        hold off
        drawnow;
    end

end

function x = SolveLU(L, U, y)

    % if A = LU
    % then x = A\y
    % can be obtained using 
    %z = L\y
    %x = U\z
    
    x = U\(L\y);

end

function v = ImageInterp(image, x, y)
% returns interpolated image intensities given real-valued image coordinates 

    ix = floor(x);
    iy = floor(y);
    fx = x-ix;
    fy = y-iy;
    
    stride = size(image,1);
    v = (1-fx).*(1-fy).*image(iy+stride.*ix) + ...
        (1-fx).*(fy).*image((iy+1)+stride.*ix) + ...
        (fx).*(1-fy).*image(iy+stride.*(ix+1)) + ...
        (fx).*(fy).*image((iy+1)+stride.*(ix+1));

end