%%%%Demo Fig.3
% 
% Author: Kaihua Zhang
% Email: cskhzhang@comp.polyu.edu.hk
% Date:  15/7/2010
% 
%--------------------------
clc;clear all;close all;

Img = 100*ones(100);
[nrow,ncol] = size(Img);
for i = 1:nrow
    for j = 1:ncol
         if (i-50).^2+(j-50).^2 < 900& (i-50).^2+(j-50).^2 > 100
            Img(i,j) = 150;
        end
    end
end
%-----------------------------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.
%-----------------------------------Level set function intialization
ic = nrow/2; jc = ncol/2; r = 40;
phi = sdf2circle(nrow,ncol, ic,jc,r);
%--------------------------------------

m = 1;
%-----------------------------------
figure(1);%imagesc(Img);colormap(gray);hold on;contour(phi,[0 0],'Color','b','LineWidth',1);
delt = .1;delt2 = 0.001; F = 0.5; nu=1;% parameters
L = 5;

%--------------------------------------
for n = 0:1600
          phi =  wextend('2D','sym',phi(L+1:end-L,L+1:end-L),L);%boundary condition
          [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;
          
          Delta_H = Delta(phi,1);
   %---------------------------RD level set evolution  
          phi = phi + delt*((nu*Kappa_g + g.*F).*absR);
          phi = phi + delt2*4*del2(phi);
   %----------------------------                   
       if mod(n,200) == 0
          figure(1);subplot(1,2,1);axis equal;
          imagesc(Img);colormap(gray);hold on;contour(phi,[0 0],'Color','b','LineWidth',1);
          iterNum=[num2str(n), ' iterations'];title(iterNum);hold off;pause(.01);  
          subplot(1,2,2);hold on;
          grid on;
          xlabel('x');ylabel('\phi(x)');  
          if n == 1600
            plot(phi(50,:),'r','LineWidth',2); title(iterNum);
          else
            plot(phi(50,:));title(iterNum);hold off;pause(.01); 
          end  
       end          
end
figure;mesh(phi);title('Final LSF \phi');
