%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% image toolbox
% MATLAB file
% 
% GeodesicActiveContour.m
% Geodesic Active Contour for image segmentation
%
%   NOTE: the current implementation is slightly different from the level
%   set implementation. there is NO reinitialization. the solution t->inf
%   is in general discontinuous and, so, it losts the sub-cell resolution
%   as in the original level set implementation.
%
% function phi = GeodesicActiveContour( A, b, tfinal )
%
%   phi' = div[ g(I) grad(phi) / |grad(phi)| ] * |grad(phi)|
%   dphi/dn |_{d Omega} = 0
%   I = observed image of size [m,n]
%   phi(t=0) = circle of radius min(m,n)/2*0.8
%
%   g(u)= 1/ (1+b |grad(I)|^2)
%
% input:  A: image file name
%         b: parameter in the edge detector
%         tfinal: artificial time
%         outfilename: output file name
% example: phi = GeodesicActiveContour( A, 0.1, 1000 );
%
% created:       07/30/2009
% author:        syleung@gmail
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function phi = GeodesicActiveContour( A, b, tfinal )

fileout = 5;    % number of intermediate output

dtout = tfinal/fileout;
time=0;
[m,n,p]=size(A);
xi=2:m+1;
yi=2:n+1;
phi=zeros(m+2,n+2);
B=zeros(m+2,n+2);
B(xi,yi)=A;

[x,y]=meshgrid(1:m,1:n);
xc=m/2;
yc=n/2;
r=min(m,n)/2*0.8;
phi(xi,yi) = sqrt((x-xc).^2+(y-yc).^2)-r;
        
% g(I)
    
C = gaussian(B,2,5);
[gn,ge,gs,gw]=g(C,b);
    
imshow(uint8(A))
hold on
contour(phi(xi,yi),[0 0],'r')

while time < tfinal
    
    % du/dn=0 on the boundary
    
    phi(1,yi)= phi(2,yi);
    phi(end,yi)= phi(end-1,yi);
    phi(xi,1)= phi(xi,2);
    phi(xi,end)= phi(xi,end-1);
    
    Bn = phi(xi,yi+1);
    Be = phi(xi+1,yi);
    Bs = phi(xi,yi-1);
    Bw = phi(xi-1,yi);
    Bc = phi(xi,yi);
    
    [Gn,Ge,Gs,Gw] = NormGrad(phi,1e-5);    
    Gn = gn./Gn;
    Ge = ge./Ge;
    Gs = gs./Gs;
    Gw = gw./Gw;    
        
    % phi^k+1 = phi^k + dt * div[ g(I) grad(phi) / |grad(phi)| ] * |grad(phi)|
    
    temp = Gn.*Bn + Ge.*Be + Gs.*Bs + Gw.*Bw ...
        - (Gn+Ge+Gs+Gw).*Bc;
    
    temp = temp.* sqrt( (Bn-Bs).^2+(Be-Bw).^2 )/2;
    
    dt = 0.45;        
    dt = min(dt,tfinal-time);    
    time=time+dt;
    [time dt tfinal]
    
    phi(xi,yi)=phi(xi,yi)+dt*temp;

    % check for intermediate output
    
    if (floor(time/dtout) ~= floor((time-dt)/dtout) )
        figure
        imshow(uint8(A))
        hold on
        contour(phi(xi,yi),[0 0],'r')
    end
    
end

phi=phi(2:end-1,2:end-1);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [gn,ge,gs,gw]=g(u,b)

% g(u)= 1/ (1+b |grad(u)|^2)

[m,n]=size(u);
xi=[2:m-1];
yi=[2:n-1];

uxp=u(xi+1,yi);
uxm=u(xi-1,yi);
uyp=u(xi,yi+1);
uym=u(xi,yi-1);
uc=u(xi,yi);

duxp=uxp-uc;
duxm=uc-uxm;
duyp=uyp-uc;
duym=uc-uym;

duxc=(duxp+duxm)/2;
duyc=(duyp+duym)/2;

gn=1./(1+b*(duxc.^2+duyp.^2));
ge=1./(1+b*(duxp.^2+duyc.^2));
gs=1./(1+b*(duxc.^2+duym.^2));
gw=1./(1+b*(duxm.^2+duyc.^2));

return

function [gn,ge,gs,gw]=NormGrad(u,eps)

% NormGrad(u)= |grad(u)|_eps

[m,n]=size(u);
xi=[2:m-1];
yi=[2:n-1];

uxp=u(xi+1,yi);
uxm=u(xi-1,yi);
uyp=u(xi,yi+1);
uym=u(xi,yi-1);
uc=u(xi,yi);

duxp=uxp-uc;
duxm=uc-uxm;
duyp=uyp-uc;
duym=uc-uym;

duxc=(duxp+duxm)/2;
duyc=(duyp+duym)/2;

gn=sqrt(duxc.^2+duyp.^2+eps^2);
ge=sqrt(duxp.^2+duyc.^2+eps^2);
gs=sqrt(duxc.^2+duym.^2+eps^2);
gw=sqrt(duxm.^2+duyc.^2+eps^2);

return