% CS 778, Final Project
% Ekrem Kocaguneli

% clear the workspace
clear;clc;
close all;

% inputImageName = 'brainImage.tif';
% inputImageName = 'long2.tif';
inputImageName = 'brain.jpg';

%%%%% follow the main steps of the algorithm
% read the image
myImage = double(imread(inputImageName));
if size(myImage,3) ~= 1
    myImage = myImage(:,:,1);
end

% create an initial active contour
myInitCont = zeros(size(myImage));
myPoint = ceil(size(myImage,1)/10);
myInitCont(myPoint:(size(myImage,1)-myPoint),myPoint:(size(myImage,1)-myPoint)) = 1;
% myInitCont=zeros(size(myImage,1),size(myImage,2));
% myInitCont(20:120,20:120) = 1;

%% initialization
myMu = 0.1;         % weight for the length term, taken as 0.2 as given in the paper
myLambda1 = 1;      % default value of 1 given in the paper is taken
myLambda2 = 1;      % default value of 1 given in the paper is taken
myEpsilon = 1;      % i took epsilon=1 as in the paper
myTimestep = 1;     % define the time-step as 0.1, as given in the paper

% initialize the level set function
myPhi0 = bwdist(myInitCont)- bwdist(1-myInitCont) + im2double(myInitCont) - 0.5;

%% the main part of the method, i.e. iterations
maxIterations = 1000;
myF = eps;

hold on;
for iterCount = 1:maxIterations
    % get indices for the inside and the outside the curve
    myInside = find(myPhi0>=0);
    myOutside = find(myPhi0<0);
    dPhi_dt = 0;
    
    % calculate the heaviside of myPhi0
    myZ = myPhi0;
    myH = zeros(size(myPhi0));
    myZ1Indices = find(myZ > myEpsilon);
    myZOtherIndices = find(myZ < myEpsilon & myZ > -myEpsilon); % absolute value case
    myH(myZ1Indices) = 1;
    for i=1:length(myZOtherIndices)
        myH(myZOtherIndices(i)) = 0.5 * ...
            (1 + myZ(myZOtherIndices(i))/myEpsilon + ... 
            (1/pi) * sin(pi * myZ(myZOtherIndices(i)) / myEpsilon) );
    end;
    
    % calculate the averages inside/outside the curve (i.e. c1 and c2)
    % below
    tmpImage = im2Double(myImage);
    myC1 = sum(sum(tmpImage .* myH))/length(myInside + eps); % take average of inside points
    myC2 = sum(sum(tmpImage .* (1 - myH)))/length(myOutside + eps); % take average of outside points
    
    % then calculate initial contour: d(phi)/dt
    dPhi_dt = -(tmpImage - myC1).^2 + (tmpImage - myC2).^2 + dPhi_dt;
    
    myF = myMu * myPhi0 ./ max(abs(myPhi0(:))) + dPhi_dt;
    myF = myF ./ (max(max(abs(myF))));
    
    
    % control if we want to stop, if not then update myF
    myPhi0New = myPhi0 + myTimestep .* myF;
    indicesToCheck = find(myPhi0New <= myTimestep);
    changeAmount = sum(abs(myPhi0(indicesToCheck) - myPhi0New(indicesToCheck)));
    if changeAmount < 60
        imagesc(myImage);colormap(gray);
        contour(myPhi0, [0 0], 'b','LineWidth',2);
        title(['Stopped Due to Too Small Change at Iteration Count: ' num2str(iterCount) ]);
        drawnow;
       break; 
    else
        myPhi0 = myPhi0New;
    end
    
    % show the image at every 10 iterations
    if mod(iterCount,10) == 0
        imagesc(myImage);colormap(gray);
        contour(myPhi0, [0 0], 'b','LineWidth',2);
        title(['Iteration Count: ' num2str(iterCount) ]);
        drawnow;
    end
end
hold off;
