%%%%Demo Fig.9
% Use flag to choose different methods
% Author: Kaihua Zhang
% Email: cskhzhang@comp.polyu.edu.hk
% Date:  15/7/2010
% Revised: 11/7/2011
% Revised: 2/8/2011
clc;clear all;close all;
%---------------
fprintf('You can choose different methods\n')
fprintf('1: RD;2: Re-initialization method 3: GDRLSE1; 4:GDRLSE2; 5:GDRLSE3 \n\n\n')
%--------------------------
flag= input('Please choose your method and then press Enter:\n 1: RD;2: Re-initialization method 3: GDRLSE1; 4:GDRLSE2; 5:GDRLSE3\n');
% flag = 1; %1: RD;2: Re-initialization method 3: GDRLSE1; 4:GDRLSE2; 5:GDRLSE3
%---------------
Img = imread('noisyImg.bmp');
Img = double(Img(:,:,1));
[nrow,ncol]=size(Img);
%---------------Intialization
%------------
switch flag
    case 1
          r = 10;c0 = 3;
          phi = c0*ones(nrow,ncol);
          phi(r:nrow-r,r:ncol-r) = -c0;
      
    case 2
          ic = nrow/2; jc = ncol/2; r = 30;
          phi = sdf2circle(nrow,ncol, ic,jc,r); 
    case 3
          r = 10;c0 = 3;
          phi = c0*ones(nrow,ncol);
          phi(r:nrow-r,r:ncol-r) = -c0;              
    case 4
          r = 10;c0 = 3;
          phi = c0*ones(nrow,ncol);
          phi(r:nrow-r,r:ncol-r) = -c0;              
    case 5
          r = 10;c0 = 3;
          phi = c0*ones(nrow,ncol);
          phi(r:nrow-r,r:ncol-r) = -c0;
    

end

imagesc(Img);colormap(gray);hold on;contour(phi,[0 0],'r','LineWidth',1);

%---------------Level set function intialization
%-----------------------------compute edge-stopping function
sigma=1.5;    % scale parameter in Gaussian kernel for smoothing. 1.5, twocells.0.45,synthetic image
G=fspecial('gaussian',15,sigma);
Img_smooth=conv2(Img,G,'same');  % smooth image by Gaussian convolution

[Ix,Iy]=gradient(Img_smooth);
f=Ix.^2+Iy.^2;
g=1./(1+f);  % edge indicator function.
%--------------------------------------------------------------------------
%--------Parameters
switch  flag
    case 1
        delt = .1;delt2 = 0.001; F =.5; nu=1;Iternum = 3000;
    case 2
        delt = .1; F = 0.5; nu = 1;Iternum = 3000;  
    case 3
        
        delt = .1; F = 0.5;mu=0.2; nu = 1;Iternum = 500;
    case 4
         delt = .1;F =.5;mu=.2; nu=1;Iternum = 5000;            
    case 5
        delt = .1; F = 0.5;mu=0.2; nu = 1;Iternum = 3000; 
end
%----------------------------------------------------
for n=1:Iternum
    phi = NeumannBoundCond(phi);
    [Kappa,absR] = Curvature_central(phi);
    %---------
    [nx,ny]=gradient(phi);
    absR = sqrt(nx.^2+ny.^2);
    absR = absR +(absR==0)*eps;
    [nxx,junk]=gradient(g.*nx./absR);  
    [junk,nyy]=gradient(g.*ny./absR);
    Kappa_g=nxx+nyy;
    %--------------------------------------
     
    switch flag
        case 1   
                     phi = phi + delt*((nu*Kappa_g + g.*F).*absR);
                     phi = phi + delt2*4*del2(phi);  
                     
                     Tit = 'RD';
        case 2   
                     phi = phi + delt*((nu*Kappa_g + g.*F).*absR); 
                     phi = reinit_SD(phi, 1, 1, 0.5, 5);
                    
                     Tit = 'Re-initialization';
        case 3       %-------------------------GDRLSE1
                     phi = phi + delt*(mu*(4*del2(phi)-Kappa)+(nu*Kappa_g + g.*F).*absR);
                     Tit = 'GDRLSE1';
        case 4       %----------------GDRLSE2 
                     r = RegularizeF(phi);
                     phi = phi + delt*(mu*r+(nu*Kappa_g + g.*F).*absR);
                     Tit = 'GDRLSE2'; 
        case 5       %----------------GDRLSE3
                     Kappa_H = Reg_Xie(phi);
                     phi = phi + delt*(mu*Kappa_H+(nu*Kappa_g + g.*F).*absR);
                     Tit = 'GDRLSE3';
                    
    end
         
   if mod(n,100)==0
       imagesc(Img);colormap(gray);hold on;contour(phi,[0 0],'b','LineWidth',1);
       iterNum=[num2str(n), ' iterations'];pause(.01);title(Tit);
   end
end
figure;mesh(phi);
title(strcat('Final LSF of',Tit));
