% This Matlab file demomstrates a variational level set method that improves the original method in Li et al's paper
%    "Level Set Evolution Without Re-initialization: A New Variational Formulation"
%    in Proceedings of CVPR'05, vol. 1, pp. 430-436.
% Author: Chunming Li, all rights reserved.
% E-mail: li_chunming@hotmail.com
% URL:  http://www.engr.uconn.edu/~cmli/

clear all;
close all;
Img=imread('noisyStar_SNR_20_to_10.bmp');
sigma=1.5;    % scale parameter in Gaussian kernel for smoothing.
G=fspecial('gaussian',15,sigma);
Img_smooth=conv2(Img,G,'same');  % smooth image by Gaussiin convolution
[Ix,Iy]=gradient(Img_smooth);
f=Ix.^2+Iy.^2;
g=1./(1+f);  % edge indicator function.

epsilon=1.5; % the papramater in the definition of smoothed Dirac function

timestep=1;  % time step
mu=0.2;  % coefficient of the internal (penalizing) energy term P(\phi)
          % Note: The product timestep*mu must be less than 0.25 for stable evolution

lambda=5; % coefficient of the weighted length term Lg(\phi)
alf=0;   % coefficient of the weighted area term Ag(\phi);
         % Note: Choose a positive(negative) alf if the initial contour is outside(inside) the object.


% define initial level set function (LSF) as -c0, c0 at points outside and inside of a region R, respectively.
[nrow, ncol]=size(Img);  
c0=2;    % The constant value used to define binary level set function as initial LSF;
         % Using initial LSF with smaller value of c0 usually speed up the evolution.

% Choose an appropriate threshold and the parameter alf. 
threshChoice=1;
if threshChoice == 1
    T=110; alf=0;   % choose the mean of all pixel intensities
elseif threshChoice == 2
    T=120; alf=-2; % too high threshold --> loss of foreground --> use negative alf to expand the contour 
else
    T=100; alf=2;  % too low threshold  --> loss of background --> use positive alf to shrink the contour 
end
% Note: the intensities of this test image are generated as Gaussian randan
% numbers with mean 120 in the foreground (the star) and 100 in the background and noise
% standard deviation 10. The above three thresholds are the mean in the
% background, foreground, and the the average of them. In practice, these
% quantities are unknown. They can be estimated by some statistical
% computation, otherwise the user need to adjust the threshold and the
% parameter alf to get desired result. Note that this method is suitable for noisy bi-modal
% images (histogram with two peaks).

BW=(Img>T);  % Define initial LSF from thresholding
initialLSF=2*c0*(0.5-BW);

u=initialLSF;

figure;
imagesc(u);
title('Initial level set function');

figure;
imagesc(Img, [0, 255]);colormap(gray);hold on;
[c,h] = contour(u,[0 0],'r');                          
title('Initial contour');                        

% start level set evolution
for n=1:60
    u=EVOLUTION_LSD(u, g ,lambda, mu, alf, epsilon, timestep, 1);      
    if mod(n,1)==0
        pause(0.01);
        imagesc(Img, [0, 255]);colormap(gray);hold on;
        [c,h] = contour(u,[0 0],'r'); 
        iterNum=[num2str(n), ' iterations'];        
        title(iterNum);
        hold off;
    end
end
imagesc(Img, [0, 255]);colormap(gray);hold on;
[c,h] = contour(u,[0 0],'r'); 
totalIterNum=[num2str(n), ' iterations'];  
title(['Final contour, ', totalIterNum]);

figure;
mesh(u);
title('Final level set function');