%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEGMENTATION.M - Detecting objects using level sets
%
% Capturing Objects in Images
%
%   Benjamin Ong, Simon Fraser University
%   Sept 14, 2001
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
close all;

% Geometry and mesh
hmesh=0.1;              % Approximate element size, course mesh

fem.geom=geomcsg({rect2(0,1,0,1)});
fem.mesh=meshinit(fem,'hmax',hmesh);
fem.mesh=meshrefine(fem);
fem.mesh=meshrefine(fem);
fem.mesh=meshrefine(fem);
fem.mesh=meshrefine(fem);

%initialize image

disp('creating image');

x = fem.mesh.p(1,:);
y = fem.mesh.p(2,:);
r = ((x-0.3).^2+(y-0.7).^2).^(0.5);
rr = ((x-0.5).^2+(y-0.2).^2).^(0.5);

[m,n] = size(x);
image = zeros(m,n);

for k = 1:n
    if ((y(k) <= 0.9) & (x(k) <= 0.9) & (y(k) >= -x(k)+1.5))
        image(k) = 255;
    end
    
    if ((r(k) <= 0.25) & (r(k) >= 0.1))
        image(k) = 255;
    end
    
    if ((rr(k) >= 0.15) &  (x(k) >= 0.3) & (x(k) <= 0.7) & (y(k) <= 0.4) & (y(k)>=0.1))
        image(k) = 255;
    end
    
end

clear k;
clear m;
clear n;
clear r;
clear rr;


image = image'; %storing it as a vector of size [m,1]


epsilon = 1e-4;
p = 2;
mu = 0.1;
nu = 0;
lambda1 = 1;
lambda2 = 1;
                                         
% Re-initialization
%to be filled in later


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FEM Structures
%   femphi - Advection of distance function phi
%   femri  - Re-initialization of distance function phi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% notice that we treat the delta term and the non linear term explicitly.  
% only the diff/grad term is treated implicitly.

disp('manually inserting initial contour');

phi = sqrt((x-0.5).^2+(y-0.5).^2)-0.3;


H = heaviside(epsilon,phi');
D = delta(epsilon,phi');

% Calculating the constant c1
temp.mesh = fem.mesh;

disp('calculating c1');

temp.sol.u = image.*H;
top_c1 = postint(temp,'u');

temp.sol.u = H;
bottom_c1 = postint(temp,'u');

disp('calculating c2');

temp.sol.u = image.*(1-H);
top_c2 = postint(temp,'u');

temp.sol.u = (1-H);
bottom_c2 = postint(temp,'u');


c1 = top_c1/bottom_c1;
c2 = top_c2/bottom_c2;

disp('initializing L to the initial contour, i.e. a circle of radius 0.3');

L = pi*0.3.^2;


disp('Setting PDE coefficients');

fem.sol.u = image;
fem.variables.image = posteval(fem,'u','tpoint','gauss2');
fem.sol.u = D';
fem.variables.D = posteval(fem,'u','tpoint','gauss2');
fem.sol = {};

fem.init = phi';

fem.variables.c1 = c1;
fem.variables.c2 = c2;
fem.variables.mu = mu;
fem.variables.lambda1 = lambda1;
fem.variables.lambda2 = lambda2;
fem.variables.p = p;
fem.variables.L = L;
fem.variables.nu = nu;


clear equ
equ.da={{{'1'}}};
equ.c = {{{'mu*p*D*L^(p-1)./sqrt(ux.^2+uy.^2)'}}};
equ.al={{{'0';'0'}}};
equ.ga={{{'0';'0'}}};
equ.be={{{'0';'0'}}};
equ.a={{{'0'}}};
equ.f = {{{'-nu - lambda1*(image-c1).^2 + lambda2*(image-c2).^2'}}};
equ.ind=1;
fem.equ=equ;

fem.dim={'u'};
fem.sdim = {'x','y'};

% Boundary conditions
fem.border=0;

% Usage
fem.usage=1;

% Problem form
fem.form='coefficient';

% Differentiation
fem.diff='off';

% Differentiation simplification
fem.simplify='on';

% Differentiation rules
fem.rules={};

% Boundary conditions
clear bnd
bnd.q={{{'0'}}};
bnd.g={{{'0'}}};
bnd.h={{{'1'}}};
bnd.r={{{'0'}}};
bnd.var={};
bnd.ind=ones(1,4);
bnd.type={'dir'};
fem.bnd=bnd;

disp('solving problem')
%fem.sol = femtime(fem,'report','on','tlist',linspace(0,0.01,0.2));

fem.sol=femtime(fem,...
	'tlist',  0:1e-6:1e-4,...
	'atol',   0.00001,...
	'rtol',   0.0001,...
	'jacobian','equ',...
	'mass',   'full',...
	'ode',    'ode15s',...
	'odeopt', struct('InitialStep',{[]},'MaxOrder',{5},'MaxStep',{[]}),...
	'out',    'sol',...
	'stop',   'on',...
	'report', 'on',...
	'context','local',...
	'sd',     'off',...
	'Epoint', 'gauss2',...
	'Nullfun','flnullorth',...
	'Tpoint', 'gauss2',...
	'Solcomp',1);


disp('displaying solution');

postmovie(fem,...
	'context','local',...
	'contdata',{'u','cont','on'},...
	'contlevels',[0 0.1 -0.1],...
	'contstyle','color',...
	'contlabel','off',...
	'contmaxmin','off',...
	'contbar','on',...
	'contmap','cool',...
	'geom',   'on',...
	'geomcol','bginv',...
	'movie',  'on',...
	'repeat', 1,...
	'fps',    1,...
	'title',  'Contour: u (u)  ',...
	'renderer','zbuffer')