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

x = 1:100;
y = 1:100;
z = 1:100;
[X,Y,Z]=meshgrid(x,y,z);
% % -------------------Generate dumbbell
offset = 20; r = 10; w = 5;
Left = sqrt((X-50+offset).^2+(Y-50).^2+(Z-50).^2)-r;
Right = sqrt((X-50-offset).^2+(Y-50).^2+(Z-50).^2)-r;
Center = max(abs(X-50)-offset,sqrt((Y-50).^2+(Z-50).^2)-w);
phi = min(Center,min(Left,Right));
%---------------------------------------
%-----------Parameters
delt = 0.1;
delt2 = 0.1;
lambda = 1;
alpha  = 0;
%--------------------------------------
fprintf('The evolution will take some minutes, please wait patiently...')
%------------------------
for n = 0:100
    strcat(num2str(n),' iterations...')
    %-------------------------------------------
    [Fx,Fy,Fz] = gradient(phi);
    deltaF = sqrt(Fx.^2+Fy.^2+Fz.^2);
    deltaF = deltaF + (deltaF==0).*eps;
    kappa = curvature_3d(phi);
    %----------------------------RD Level set evolution
    phi = phi + delt*(lambda*kappa+alpha).*deltaF;
    phi = phi + delt2*4*del2(phi);
    %----------------------------
    if mod(n,20)==0
    p = patch(isosurface(X,Y,Z,phi,0));
    isonormals(X,Y,Z,phi,p);
    set(p,'FaceColor','red','EdgeColor','none');
    daspect([1 1 1]);
    view(3); 
    camlight ;axis tight;
    lighting gouraud  
    iterNum=[num2str(n), ' iterations'];title(iterNum);
    pause(0.1);
    hold off;
    clf('reset'); 
    end
end
    p = patch(isosurface(X,Y,Z,phi,0));
    isonormals(X,Y,Z,phi,p);
    set(p,'FaceColor','red','EdgeColor','none');
    daspect([1 1 1]);
    view(3); 
    camlight ;axis tight;
    lighting gouraud  
    iterNum=[num2str(n), ' iterations'];title(iterNum);
    pause(0.1);
    hold off;