% Run with for instance:
% [phi1,phi2,psi,segm,binwat,F] = segmsurftopo(im,minima,prm);                                    
%
% Use in combination with GETIM in the folder under
%
function [phi1,phi2,psi,segm,F] = segmsurftopo(varargin)

tic

global M N methodTdf

im = varargin{1};
minima = varargin{2};
prmin = varargin{3};
for i = 4 : 2 : nargin
    varname = varargin{i};
    varhere = varargin{i+1};
    switch(varname)
        case('topo')
            c = varhere;
        otherwise
            error('Wrong input')
    end;
end;



% solver
% solver = 1 : Chan-Vese
% solver = 2 : PCLS (Standard from polynomials)
% solver = 3 : Bin
% solver = 4 : PCLS2 (time step 0.05-0.1,gamma=0.01,epsPCLS=1e-01,lambda=0.05)
% solver = 5 : PCLS3

% default values
prm.gamma = 0.02;
prm.epsPCLS = 0.02;
prm.lambda = 0.3;
prm.vis = Inf;
prm.dt = 0.002;
prm.maxiter = 500;
% the method fr regularsation
% methodreg = 1 : on phi
% methodreg = 2 : on psi
prm.methodreg = 2;
prm.solver = 3;

% assign the input variables
if ~isempty(prmin)
    names = fieldnames(prmin);
    numfields = size(names,1);
    for i = 1 : numfields
        namehere = names{i};
        valfield = eval(['prmin.' namehere]);
        prm = setfield(prm,namehere,valfield);
    end;
end;
% to make the old stuff work.....
methodreg = prm.methodreg;
dt = prm.dt;
maxiter = prm.maxiter;
solver = prm.solver;

if solver == 1
    namesolver = 'ChanVese';
elseif solver == 2
    namesolver = 'PCLS polynomial';
elseif solver == 3
    namesolver = 'Binary';
elseif solver == 4
    namesolver = 'PCLS, new';
elseif solver == 6
    namesolver = 'PCLS, two phases, then four';
else
    error('Wrong solver')
end;

% lambda:
% reguarization term in Chan-Vese and binary
% liten ved method == 1, ellers 0.1, not too small eps in psi definition,
% then too singular

[M N] = size(im);

% get parameters
[gamma,epsPCLS,lambda,vis,sigma] = getprm(solver,prm);

im = scale(im);


% stop iteration if change drops below
thCh = 0.00000;

% method for TDF
% METHODTDF = 1 : IFT
% METHODTDF = 2 : Flooding
% does this work???
methodTdf = 1;


%
% Get image, minima and S
%

% find boundary and color image to calculate S
[partition,boundim] = distim(minima,8,1);
colorim = findCol(boundim);
    
p = 1;
imred = imresize(im,p,'bicubic');
colorimred = imresize(colorim,p,'nearest');


% find topographical distance function
C = tdf(imred,colorimred);

for i = 1 : 4
    c(:,:,i) = imresize(C(:,:,i),[M N],'bicubic'); 
end;

% check size of c
[m n o] = size(c);
for i = o+1 : 4
    c(:,:,i) = ones(M,N);
end;


% disp(sprintf('Name of experiment: %s',name));

[M N] = size(im);

%
% Initilaize phi
%
[phi1 phi2] = iniphi(solver,colorim,c);


if maxiter == 0
    % Basis functions to plot
    psi = getpsi(phi1,phi2,solver);

    % Make segmentation
    segm = getsegm(psi);

    % make watershed lines
    binwat = getbinwat(segm);       
    
    F = 0;
    return;
end;

disp(sprintf('Using solver %s',namesolver));

% main iteration
niter = 0;F = 0;
while 1

    niter = niter + 1;
    
   
    % for F
    phiOld1 = phi1;
    phiOld2 = phi2;
    
    % solve PDE    
    if solver == 1
        [phi1,phi2,F(niter)]= solveFourPhasesCVFourCol(c,phi1,phi2,lambda,dt,methodreg);
    elseif solver == 2
        [phi1,F(niter)]= solveFourPhasesPCLSFourCol(c,phi1,colorIm,boundIm,dt);
        phi2 = phi1;
    elseif solver == 3
        [phi1,phi2,F(niter)]= solveFourPhasesBinFourCol(c,phi1,phi2,dt,methodreg,sigma,lambda);    
    elseif solver == 4
        [phi1,F(niter)]= solveFourPhasesPCLS2FourCol(c,phi1,dt,lambda,methodreg,epsPCLS,gamma);
        phi2 = phi1;                                
    elseif solver == 5
        [phi1,F(niter)]= solveFourPhasesPCLS3(c,phi1,dt);
         phi2 = phi1;
    elseif solver == 6
        if niter < 150
            [phi1,F(niter)]= solveFourPhasesPCLS5FourCol(c,phi1,dt,lambda,methodreg,epsPCLS,gamma);
        else
            [phi1,F(niter)]= solveFourPhasesPCLS2FourCol(c,phi1,dt,lambda,methodreg,epsPCLS,gamma);
        end;
        phi2 = phi1;                                
    end;
    

    % find max change    
    ch1 = max(max(abs (phiOld1 - phi1)));
    ch2 = max(max(abs (phiOld2 - phi2)));
    ch = max(ch1,ch2);
   
    
%     % if we proceed towards end of iterations then switch on constraint
%     if niter > 0.7*maxiter | ch < 2*thCh
%         gamma = gammaFinal;
%         epsPCLS = epsPCLSFinal;
%     end;
    
    disp(sprintf('Number of iterations : %i',niter));
%     disp(sprintf('Maximum change : %6.5f',ch));
%     disp(sprintf('Gamma : %i',gamma));
%     disp(sprintf('epsPCLS : %6.5f',epsPCLS));
%     disp(sprintf('sigma : %6.5f',sigma));
%     disp(sprintf('lambda : %6.5f',lambda));
%     disp(sprintf('dt : %6.5f',dt));

    % test if phi changes
    if ch < thCh
        break;
    end;
    
    % test if max number of iterations
    if niter == maxiter
        break;
    end;
    
    if round(niter/vis) == niter/vis && vis < inf
        show(phi1,30);colorbar;title('phi1');drawnow;
        show(phi2,31);colorbar;title('phi2');drawnow;
    end;
    
end;

%
% Basis functions to plot
%
psi = getpsi(phi1,phi2,solver);

% showall(psi)

%
% Make segmentation
%
segm = getsegm(psi);


if vis < inf
    figure(20);
    subplot(3,2,1);imagesc(psi(:,:,1));colormap(gray);colorbar;axis off;title('psi1');
    subplot(3,2,2);imagesc(psi(:,:,2));colormap(gray);colorbar;axis off;title('psi2');
    subplot(3,2,3);imagesc(psi(:,:,3));colormap(gray);colorbar;axis off;title('psi3');
    subplot(3,2,4);imagesc(psi(:,:,4));colormap(gray);colorbar;axis off;title('psi4');
    subplot(3,2,5);imagesc(im);colormap(gray);colorbar;title('im');axis off        
    subplot(3,2,6);imagesc(segm);colormap(gray);colorbar;title('segm');axis off
end;

%
% Make watershed lines
%
binwat = getbinwat(segm);

% show(im,1);colorbar
% show(phi1,6);colorbar;title('phi1')
% show(phi2,7);colorbar;title('phi2')
% show(segm,8);colorbar
% show(binWat,9);colorbar;title('Watershed levelset')
% show(binWatStand,10);colorbar;title('Standard watershed, watersheds')
% show(watStand,11);colorbar;title('Standard watershed, regions')
% figure(12);plot(F);

toc

%----------------------------------------------

function [binwat] = getbinwat(segm)

global M N;

binwat = zeros(M,N);
n = M*N;
for i = 1 : max(segm(:))
    regHere = eq(segm,i);
    ind = find(regHere);
    
    if size(find(regHere),1) < 100 
        continue;
    end
    
    perim = bwperim(regHere);
    binwat(perim) = 1;

end;
binwat = bwmorph(binwat,'thin','Inf');
binwat = constborder(binwat,1,0);

%------------------------------------------------


function[phi1 phi2] = iniphi(solver,colorim,c)

global M N O

% initialize level set functions
if solver == 1
    phi1 = zeros(M,N);  
    phi2 = zeros(M,N);  
elseif solver == 2
    phi1 = 2*ones(M,N);
    phi2 = phi1;
elseif solver == 3
    phi1 = zeros(M,N);
    phi2 = zeros(M,N);
elseif solver == 4 || solver == 6
    phi1 = zeros(M,N);
%    phi1 = round(rand(M,N)) + 2;
%    phi1 = imfilter(colorIm,fspecial('gaussian',5,1),'replicate');
    phi1 = 2.5*ones(M,N);
%     phi1 = 4*ones(M,N);

    
%     for i = 1 : 4
%         chere = c(:,:,i);
%         show(chere,1)
%         if isempty(find(ne(chere,1)))
%             continue;
%         end;
%         for j = 1 : 4
%             colhere = eq(colorim,j);
%             if isempty(find(colhere))
%                 val(j) = inf;
%             else
%                 val(j) = mean(chere(colhere));
%             end;
% %             val(j)
% %         show(colhere,3)
% %         pause
%         end;
%         [minval,ind] = min(val);
% %         ind
% 
%         phi1(eq(colorim,ind)) = i;
% %         show(phi1,2)
% %         pause
%     end;
    phi1(colorim <= 2) = 2;
    phi1(colorim >= 3) = 3;
%     phi1 = colorim;
    phi2 = phi1;
elseif solver == 5
    phi1 = 4*ones(M,N);
    phi2 = phi1;
end;

%-------------------------------------------

function [psi] = getpsi(phi1,phi2,solver);

if solver == 1
    eps = 0.001;
    H1 = ge(phi1,0);
    H2 = ge(phi2,0);
%     H1 = heaviside(phi1,eps);
%     H2 = heaviside(phi2,eps);
    
    psi(:,:,1) = H1.*H2;
    psi(:,:,2) = H1.*(1-H2);
    psi(:,:,3) = (1-H1).*H2;
    psi(:,:,4) = (1-H1).*(1-H2);
    psi = round(psi);
elseif solver == 2
%    phi1 = round(phi1);
    phi2 = phi1;
    phi = phi1;
    % basis functions exact
    alpha(1) = (1-2)*(1-3)*(1-4);
    psi(:,:,1) = (1/alpha(1))*(phi-2).*(phi-3).*(phi-4);
    
    alpha(2) = (2-1)*(2-3)*(2-4);
    psi(:,:,2) = (1/alpha(2))*(phi-1).*(phi-3).*(phi-4);
    
    alpha(3) = (3-1)*(3-2)*(3-4);
    psi(:,:,3) = (1/alpha(3))*(phi-1).*(phi-2).*(phi-4);
    
    alpha(4) = (4-1)*(4-2)*(4-3);
    psi(:,:,4) = (1/alpha(4))*(phi-1).*(phi-2).*(phi-3);    

   psi = gt(psi,0.1);
    

elseif solver == 3
    eps = 10e-6;
    absPhi1 = abs(phi1); absPhi2 = abs(phi2);
    p1 = phi1./(absPhi1 + eps); p2 = phi2./(absPhi2 + eps);
    psi(:,:,1) =   (1/4)*(1+p1).*(1+p2);
    psi(:,:,2) =   (1/4)*(1+p1).*(1-p2);
    psi(:,:,3) =   (1/4)*(1-p1).*(1+p2);
    psi(:,:,4) =   (1/4)*(1-p1).*(1-p2);
    psi = gt(psi,0.05);
elseif solver == 4 || solver == 6
    eps = 10e-6;

%     phi1 = round(phi1);
    phi2 = phi1;
    for i = 1 : 4
        psi(:,:,i) = 0.5*((phi1-i+0.5)./sqrt((phi1-i+0.5).^2 + eps) ...
            - (phi1-i-0.5)./sqrt((phi1-i-0.5).^2 + eps));
    end;
    psi = gt(psi,0.3);
%     show(phi1,1)
%     show(psi(:,:,1),2)
%     pause
elseif solver == 5
    eps = 10e-6;
    numObj = size(c,3);
    for i = 1 : numObj
        psi(:,:,i) = 0.5*((phi1-i+0.5)./sqrt((phi1-i+0.5).^2 + eps) ...
            - (phi1-i-0.5)./sqrt((phi1-i-0.5).^2 + eps));
    end;
    psi = gt(psi,0.2);
end;

%-------------------------------------------------------

function [segm] = getsegm(psi)

global M N

segm = zeros(M,N);
teller = 0;

for i = 1 : 4
    psiHere = psi(:,:,i);
%    ind = find(psiHere > 0.5);
%    psiHere(ind) = 1;
%     show(psiHere,3);colorbar    
    [faser L] = bwlabeln(psiHere);
    
    for j = 1 : L
        teller = teller + 1;
        regHere = eq(faser,j);
        segm(regHere) = teller;            
    end;

end;

%----------------------------------------------------


% polynomial
function [phi,F]= solveFourPhasesPCLSFourCol(c,phi,colorIm,boundIm,dt);

global method gamma


[M N] = size(phi);

% basis functions
alpha(1) = (1-2)*(1-3)*(1-4);
psi(:,:,1) = (1/alpha(1))*(phi-2).*(phi-3).*(phi-4);

alpha(2) = (2-1)*(2-3)*(2-4);
psi(:,:,2) = (1/alpha(2))*(phi-1).*(phi-3).*(phi-4);

alpha(3) = (3-1)*(3-2)*(3-4);
psi(:,:,3) = (1/alpha(3))*(phi-1).*(phi-2).*(phi-4);

alpha(4) = (4-1)*(4-2)*(4-3);
psi(:,:,4) = (1/alpha(4))*(phi-1).*(phi-2).*(phi-3);


% derivatives
dPsi1_dPhi = (1/alpha(1))*((phi-3).*(phi-4)+(phi-2).*(phi-4)+(phi-2).*(phi-3));
dPsi2_dPhi = (1/alpha(2))*((phi-3).*(phi-4)+(phi-1).*(phi-4)+(phi-1).*(phi-3));
dPsi3_dPhi = (1/alpha(3))*((phi-2).*(phi-4)+(phi-1).*(phi-4)+(phi-1).*(phi-2));
dPsi4_dPhi = (1/alpha(4))*((phi-2).*(phi-3)+(phi-1).*(phi-3)+(phi-1).*(phi-2));


% constraint
K = (phi-1).*(phi-2).*(phi-3).*(phi-4);
dK_dPhi = (phi-2).*(phi-3).*(phi-4)+(phi-1).*(phi-3).*(phi-4)+ ...
    (phi-1).*(phi-2).*(phi-4)+(phi-1).*(phi-2).*(phi-3);

% update Euler Lagrange equation
update = (c(:,:,1) .* dPsi1_dPhi + c(:,:,2) .* dPsi2_dPhi + ...
          c(:,:,3) .* dPsi3_dPhi + c(:,:,4) .* dPsi4_dPhi) + (2/gamma)*K.*dK_dPhi;

phi = phi - dt*update;

F = c(:,:,1).*psi(:,:,1) + c(:,:,2).*psi(:,:,2) ...
  + c(:,:,3).*psi(:,:,3) + c(:,:,4).*psi(:,:,4);
F = sum(F(:));

% make latex file

%---------------------------------------------------

% my own ide to skip 4 color and use many more phases
function [phi,F]= solveFourPhasesPCLS3(c,phi,dt);

global method gamma

[M N] = size(phi);

eps = 1e-02;

% showall(c,'colorbar');

numObj = size(c,3);

% basis functions
for i = 1 : numObj
    psi(:,:,i) = ((phi-i+0.5)./sqrt((phi-i+0.5).^2 + eps) ...
                - (phi-i-0.5)./sqrt((phi-i-0.5).^2 + eps));
end;

% derivatives
for i = 1 : numObj
    % teller 1
    t1 = phi-i+0.5;
    % nevner 1
    n1 = sqrt((phi-i+0.5).^2 + eps);
    % teller 2
    t2 = phi-i-0.5;
    % nevner 2
    n2 = sqrt((phi-i-0.5).^2 + eps);    
    % derivative
    dPsi_dPhi(:,:,i) = ((n1.^2 - t1.^2)./(n1.^3) - (n2.^2 - t2.^2)./(n2.^3));
end;

% update Euler Lagrange equation
update = zeros(M,N);
for i = 1 : numObj
    update = update + c(:,:,i) .* dPsi_dPhi(:,:,i);
end;
% update = update + (2/gamma)*K.*dK_dPhi;
% u = max(update(:))
% show(update,3);colorbar

% update PDE
phi = phi - dt*update;

F = zeros(M,N);
for i = 1 : numObj
    F = F + c(:,:,i).*psi(:,:,i);
end;
F = sum(F(:));

%----------------------------------------------------


% the new method of Tai
function [phi,F]= solveFourPhasesPCLS2FourCol(c,phi,dt,lambda,methodreg,epsPCLS,gamma)
                 
global M N niter


[M N] = size(phi);

% basis functions
for i = 1 : 4
    psi(:,:,i) = 0.5*((phi-i+0.5)./sqrt((phi-i+0.5).^2 + epsPCLS) - ...
                  (phi-i-0.5)./sqrt((phi-i-0.5).^2 + epsPCLS));
end;


% derivatives
for i = 1 : 4
    % nevner 1
    n1 = sqrt((phi-i+0.5).^2 + epsPCLS);
    % nevner 2
    n2 = sqrt((phi-i-0.5).^2 + epsPCLS);    
    % derivative
    dPsi_dPhi(:,:,i) = 0.5*(epsPCLS ./ (n1.^3) - epsPCLS ./ (n2.^3));
    
%    dPsi_dPhi(:,:,i) = ((n1.^2 - t1.^2)./(n1.^3) - (n2.^2 - t2.^2)./(n2.^3));
end;


% constraint
K = (phi-1).*(phi-2).*(phi-3).*(phi-4);
dK_dPhi = (phi-2).*(phi-3).*(phi-4)+(phi-1).*(phi-3).*(phi-4)+ ...
    (phi-1).*(phi-2).*(phi-4)+(phi-1).*(phi-2).*(phi-3);

% update Euler Lagrange equation data term
update = zeros(M,N);
for i = 1 : 4
    updatehere = c(:,:,i) .* dPsi_dPhi(:,:,i);
    update = update + updatehere;
    if niter > Inf
    show(updatehere,1);colorbar
    show(c(:,:,i),2);colorbar
    show(psi(:,:,i),3);colorbar
    show(dPsi_dPhi(:,:,i),4);colorbar
    show(phi,5);colorbar
    show((2*gamma)*K.*dK_dPhi,6);colorbar
    pause
    end;
end;


% show(update,2);colorbar

% regularizer
if methodreg == 1
    update = update - lambda*curv(phi);
elseif methodreg == 2
    for i = 1 : 4
        update = update - lambda*curv(psi(:,:,i)) .* dPsi_dPhi(:,:,i);
    end;
end;
    

% show(lambda*curv(psi(:,:,i)) .* dPsi_dPhi(:,:,i),3);colorbar
% pause
update = update + (2*gamma)*K.*dK_dPhi;
% u = max(update(:))
% show(update,3);colorbar

phi = phi - dt*update;

F = c(:,:,1).*psi(:,:,1) + c(:,:,2).*psi(:,:,2) ...
  + c(:,:,3).*psi(:,:,3) + c(:,:,4).*psi(:,:,4);
F = sum(F(:));

%------------------------------------------------------------------
% the new method of Tai, first converging towards 2 and 3, then converging
% towards 1,2,3,4 since the pixels can not cross another phase
function [phi,F]= solveFourPhasesPCLS5FourCol(c,phi,dt,lambda,methodreg,epsPCLS,gamma)
                 
global M N niter


[M N] = size(phi);

% basis functions
for i = 2 : 3
    psi(:,:,i) = 0.5*((phi-i+0.5)./sqrt((phi-i+0.5).^2 + epsPCLS) - ...
                  (phi-i-0.5)./sqrt((phi-i-0.5).^2 + epsPCLS));
end;


% derivatives
for i = 2 : 3
    % nevner 1
    n1 = sqrt((phi-i+0.5).^2 + epsPCLS);
    % nevner 2
    n2 = sqrt((phi-i-0.5).^2 + epsPCLS);    
    % derivative
    dPsi_dPhi(:,:,i) = 0.5*(epsPCLS ./ (n1.^3) - epsPCLS ./ (n2.^3));
    
%    dPsi_dPhi(:,:,i) = ((n1.^2 - t1.^2)./(n1.^3) - (n2.^2 - t2.^2)./(n2.^3));
end;


% constraint
K = (phi-2).*(phi-3);
dK_dPhi = (phi-3) + (phi-2);

% update Euler Lagrange equation data term
update = zeros(M,N);
% show(c(:,:,1),1)
% show(c(:,:,2),2)
% show(c(:,:,3),3)
% show(c(:,:,4),4)
c(:,:,2) = min(c(:,:,1:2),[],3);
c(:,:,3) = min(c(:,:,3:4),[],3);
% show(c(:,:,2),5)
% show(c(:,:,3),6)
% pause
for i = 2 : 3
    updatehere = c(:,:,i) .* dPsi_dPhi(:,:,i);
    update = update + updatehere;
    if niter > Inf
        i
    show(updatehere,1);colorbar
    show(c(:,:,i),2);colorbar
    show(psi(:,:,i),3);colorbar
    show(dPsi_dPhi(:,:,i),4);colorbar
    show(phi,5);colorbar
    show((2*gamma)*K.*dK_dPhi,6);colorbar
    show(c(:,:,2),7)
    show(c(:,:,3),8)
    pause
    end;
end;


% show(update,2);colorbar

% regularizer
if methodreg == 1
    update = update - lambda*curv(phi);
elseif methodreg == 2
    for i = 2 : 3
        update = update - lambda*curv(psi(:,:,i)) .* dPsi_dPhi(:,:,i);
    end;
end;
    

% show(lambda*curv(psi(:,:,i)) .* dPsi_dPhi(:,:,i),3);colorbar
% pause
update = update + (2*gamma)*K.*dK_dPhi;

phi = phi - dt*update;

F = c(:,:,2).*psi(:,:,2) + c(:,:,3).*psi(:,:,3);
F = sum(F(:));

%---------------------------------------------------------------

% the new method of Tai, with basis functions going up at the ends,
% not working
function [phi,F]= solveFourPhasesPCLS4FourCol(c,phi,colorIm,boundIm,dt);

global method gamma lambda epsPCLS


[M N] = size(phi);
numcol = 4;

% basis functions
psi(:,:,1) = -0.5*(phi-1-0.5)./sqrt((phi-1-0.5).^2 + epsPCLS) + 0.5;
for i = 2 : numcol-1
    psi(:,:,i) = 0.5*((phi-i+0.5)./sqrt((phi-i+0.5).^2 + epsPCLS) - (phi-i-0.5)./sqrt((phi-i-0.5).^2 + epsPCLS));
end;
psi(:,:,numcol) = 0.5*(phi-numcol+0.5)./sqrt((phi-numcol+0.5).^2 + epsPCLS) + 0.5;


% derivatives
n2 = -sqrt((phi-1+0.5).^2 + epsPCLS);    
dPsi_dPhi(:,:,1) = 0.5*epsPCLS ./ (n2.^3);
for i = 2 : numcol-1
    % nevner 1
    n1 = sqrt((phi-i+0.5).^2 + epsPCLS);
    % nevner 2
    n2 = sqrt((phi-i-0.5).^2 + epsPCLS);    
    % derivative
    dPsi_dPhi(:,:,i) = 0.5*(epsPCLS ./ (n1.^3) - epsPCLS ./ (n2.^3));
end;
% nevner 1
n1 = sqrt((phi-numcol+0.5).^2 + epsPCLS);
dPsi_dPhi(:,:,numcol) = 0.5*epsPCLS ./ (n1.^3);


% constraint
K = (phi-1).*(phi-2).*(phi-3).*(phi-4);
dK_dPhi = (phi-2).*(phi-3).*(phi-4)+(phi-1).*(phi-3).*(phi-4)+ ...
    (phi-1).*(phi-2).*(phi-4)+(phi-1).*(phi-2).*(phi-3);

% update Euler Lagrange equation data term
update = zeros(M,N);
for i = 1 : 4
    update = update + c(:,:,i) .* dPsi_dPhi(:,:,i);
end;

% regularizer
if method == 1
    update = update - lambda*curv(phi);
elseif method == 2
    for i = 1 : 4
        update = update - lambda*curv(psi(:,:,i)) .* dPsi_dPhi(:,:,i);
    end;
end;
    

update = update + (2*gamma)*K.*dK_dPhi;
% u = max(update(:))
% show(update,3);colorbar

phi = phi - dt*update;

F = c(:,:,1).*psi(:,:,1) + c(:,:,2).*psi(:,:,2) ...
  + c(:,:,3).*psi(:,:,3) + c(:,:,4).*psi(:,:,4);
F = sum(F(:));

% make latex file

%---------------------------------------------------


function [phi1,phi2,F]= solveFourPhasesBinFourCol(c,phi1,phi2,dt,methodreg,sigma,lambda)


[M N] = size(phi1);

% basis functions
% eps = 0.01;
eps = 10e-3;
% absPhi1 = abs(phi1); absPhi2 = abs(phi2);
% p1 = phi1./(absPhi1 + eps); 
% p2 = phi2./(absPhi2 + eps);
p1 =  phi1./sqrt(phi1.^2 + eps); 
p2 =  phi2./sqrt(phi2.^2 + eps); 
psi1 =   (1/4)*(1+p1).*(1+p2);
psi2 =   (1/4)*(1+p1).*(1-p2);
psi3 =   (1/4)*(1-p1).*(1+p2);
psi4 =   (1/4)*(1-p1).*(1-p2);

%
% Derivatives
%

c1 = c(:,:,1);
c2 = c(:,:,2);
c3 = c(:,:,3);
c4 = c(:,:,4);

% gamma = 1000;
gamma = 0;
diamlow = 30;
numpixlow = pi*(diamlow/2)^2;

A1 = sum(psi1(:));
A2 = sum(psi2(:));
A3 = sum(psi3(:));
A4 = sum(psi4(:));
% numpixlow
% A1
% A2
% A3
% A4
% show(psi1,1)
% show(psi2,2)
% show(psi3,3)
% show(psi4,4)

epshere = 1000*numpixlow;
% dfdpsi1 = c1;
% dfdpsi2 = c2;
% dfdpsi3 = c3;
% dfdpsi4 = c4;
ddim1 = -gamma*0.5*epshere./((A1 - numpixlow).^2 + epshere).^(3/2); 
ddim2 = -gamma*0.5*epshere./((A2 - numpixlow).^2 + epshere).^(3/2);
ddim3 = -gamma*0.5*epshere./((A3 - numpixlow).^2 + epshere).^(3/2); 
ddim4 = -gamma*0.5*epshere./((A4 - numpixlow).^2 + epshere).^(3/2); 
% ddim1
% ddim2
% ddim3 
% ddim4
% pause
dfdpsi1 = c1 - ddim1;
dfdpsi2 = c2 - ddim2;
dfdpsi3 = c3 - ddim3;
dfdpsi4 = c4 - ddim4;
%   pause
% derivative of the factor inside brackets
d1 = eps./(phi1.^2 + eps).^(3/2);
d2 = eps./(phi2.^2 + eps).^(3/2);

dpsi1_dphi1 =  (1/4)*d1.*(1+p2);
dpsi1_dphi2 =  (1/4)*d2.*(1+p1);

dpsi2_dphi1 =  (1/4)*d1.*(1-p2);
dpsi2_dphi2 = -(1/4)*d2.*(1+p1);

dpsi3_dphi1 = -(1/4)*d1.*(1+p2);
dpsi3_dphi2 =  (1/4)*d2.*(1-p1);

dpsi4_dphi1 =  -(1/4)*d1.*(1-p2);
dpsi4_dphi2 =  -(1/4)*d2.*(1-p1);

% % an often used factor
% d1 = (1 - phi1.^2)./((phi1.^2 + eps).^(1/2));
% d2 = (1 - phi2.^2)./((phi2.^2 + eps).^(1/2));
% dpsi1_dphi1 = (1/4)*d1.*(1+p2); 
% dpsi1_dphi2 = (1/4)*(1+p1).*d2; 
% dpsi2_dphi1 = (1/4)*d1.*(1-p2); 
% dpsi2_dphi2 = -(1/4)*(1+p1).*d2; 
% dpsi3_dphi1 = -(1/4)*d1.*(1+p2); 
% dpsi3_dphi2 = (1/4)*(1-p1).*d2; 
% dpsi4_dphi1 = -(1/4)*d1.*(1-p2); 
% dpsi4_dphi2 = -(1/4)*(1-p1).*d2; 

%
% Regularization terms
%

if methodreg == 1
    reg1 = curv(phi1);
    reg2 = curv(phi2);
elseif methodreg == 2
    curv1 = curv(psi1);
    curv2 = curv(psi2);
    curv3 = curv(psi3);
    curv4 = curv(psi4);
    
    % for phi1
    reg1 = curv1.*dpsi1_dphi1 + curv2.*dpsi2_dphi1 + ...
        curv3.*dpsi3_dphi1 + curv4.*dpsi4_dphi1;
    
    % for phi2
    reg2 = curv1.*dpsi1_dphi2 + curv2.*dpsi2_dphi2 + ...
        curv3.*dpsi3_dphi2 + curv4.*dpsi4_dphi2;
end;

% constraint
dKdphi1 = 4*(phi1.^2-1).*phi1;
dKdphi2 = 4*(phi2.^2-1).*phi2; 

% the update terms, derivatives and chain rule
chere = dfdpsi1.*dpsi1_dphi1 + dfdpsi2.*dpsi2_dphi1 + ...
           dfdpsi3.*dpsi3_dphi1 + dfdpsi4.*dpsi4_dphi1;
update1 = chere + sigma*dKdphi1 - lambda*reg1;

chere = dfdpsi1.*dpsi1_dphi2 + dfdpsi2.*dpsi2_dphi2 + ...
        dfdpsi3.*dpsi3_dphi2 + dfdpsi4.*dpsi4_dphi2; 
update2 = chere + sigma*dKdphi2 - lambda*reg2;

% % sigma = 1;
% % update Euler Lagrange equation
% p1 = 2;
% p2 = 2;
% p3 = 2;
% p4 = 2;
% 
% chere = c1 .* (c1-p1) .* dPsi1_dPhi1 + c2 .* (c2-p2) .* dPsi2_dPhi1 + ...
%         c3 .* (c3-p3) .* dPsi3_dPhi1 + c4 .* (c4-p4) .* dPsi4_dPhi1;
% update1 =  chere + sigma*dKdPhi1 - lambda*reg1;
% 
% chere = c1 .* (c1-p1) .* dPsi1_dPhi2 + c2 .* (c2-p2) .* dPsi2_dPhi2 + ...
%         c3 .* (c3-p3) .* dPsi3_dPhi2 + c4 .* (c4-p4) .* dPsi4_dPhi2; 
% update2 = chere + sigma*dKdPhi2 - lambda*reg2;

phi1 = phi1 - dt*update1;
phi2 = phi2 - dt*update2;
% phi1 = phi1 - dt*update1;
% phi2 = phi2 - dt*update2;

F = c1.*psi1 + c2.*psi2 ...
  + c3.*psi3 + c4.*psi4;
F = sum(F(:));

% show(update1,1);colorbar
% show(update2,2);colorbar
% pause

%----------------------------------------------------


function [phi1,phi2,F]= solveFourPhasesCVFourCol(c,phi1,phi2,lambda,dt,methodreg);


[M N] = size(phi1);

% heaviside and delta
% first gt, second ge
H1 = gt(phi1,0);
H2 = gt(phi2,0);
eps = 0.01;
% H1 = heaviside(phi1,eps);
% H2 = heaviside(phi2,eps);
d1 = delta(phi1,eps);
d2 = delta(phi2,eps);

% basis functions
% % H1 = gt(phi1,0);
% H2 = gt(phi2,0);
psi1 = H1.*H2;
psi2 = H1.*(1-H2);
psi3 = (1-H1).*H2;
psi4 = (1-H1).*(1-H2);


% glatt phi som vanlig
if methodreg == 1

    % phi1
    update1 = c(:,:,1).*d1.*H2 + c(:,:,2).*d1.*(1-H2) - ...
              c(:,:,3).*d1.*H2 - c(:,:,4).*d1.*(1-H2);
    phi1 = phi1 - dt*(update1 - lambda * curv(phi1));
    
    % phi2
    update2 = c(:,:,1).*H1.*d2 - c(:,:,2).*H1.*d2 + ...
              c(:,:,3).*(1-H1).*d2 - c(:,:,4).*(1-H1).*d2;
    phi2 = phi2 - dt*(update2 - lambda * curv(phi2));

% glatt psi
elseif methodreg == 2

    curv1 = curv(psi1);
    curv2 = curv(psi2);
    curv3 = curv(psi3);
    curv4 = curv(psi4);

    % phi1
    update1 = c(:,:,1).*d1.*H2 + c(:,:,2).*d1.*(1-H2) - ...
              c(:,:,3).*d1.*H2 - c(:,:,4).*d1.*(1-H2);
          
    reg = d1.*H2.*curv1 + d1.*(1-H2).*curv2 ...
          - d1.*H2.*curv3 - d1.*(1-H2).*curv4;
    
    phi1 = phi1 - dt*(update1 - lambda * reg);
    
    % phi2
    update2 = c(:,:,1).*H1.*d2 - c(:,:,2).*H1.*d2 + ...
              c(:,:,3).*(1-H1).*d2 - c(:,:,4).*(1-H1).*d2;
    
    reg = H1.*d2.*curv1 - H1.*d2.*curv2 ...
          + (1-H1).*d2.*curv3 - (1-H1).*d2.*curv4;
    
    phi2 = phi2 - dt*(update2 - lambda * reg);
    
elseif methodreg == 3
    eps = 0.1;
    B1 = psi1./sqrt(psi1.^2 + eps);
    B2 = psi2./sqrt(psi2.^2 + eps);
    B3 = psi3./sqrt(psi3.^2 + eps);
    B4 = psi4./sqrt(psi4.^2 + eps);
    
    % phi1
    update1 = c(:,:,1).*B1.*d1.*H2 + c(:,:,2).*B2.*d1.*(1-H2) - ...
              c(:,:,3).*B3.*d1.*H2 - c(:,:,4).*B4.*d1.*(1-H2);
    phi1 = phi1 - dt*(update1 - lambda*curv(phi1));
    
    
    % phi2
    update2 = c(:,:,1).*B1.*H1.*d2 - c(:,:,2).*B2.*H1.*d2 + ...
              c(:,:,3).*B3.*(1-H1).*d2 - c(:,:,4).*B4.*(1-H1).*d2;
    phi2 = phi2 - dt*(update2 - lambda*curv(phi2));

end;
% smoothnesse term
[fx fy] = gradient(phi1);
absGrad1 = sqrt(fx.*fx + fy.*fy);
[fx fy] = gradient(phi2);
absGrad2 = sqrt(fx.*fx + fy.*fy);

F = c(:,:,1).*psi1 + c(:,:,2).*psi2 + c(:,:,3).*psi3 + c(:,:,4).*psi4 + ...
    lambda*(absGrad1 + absGrad2);
F = sum(F(:));


%-------------------------------------------------------

function [c] = mktdf(im,minima,colorIm,scTdf,name)

[M N O] = size(im);

imReg = filterRidges(im,50);
imReg = imresize(imReg,scTdf,'bicubic');
minimaReg = imresize(minima,scTdf,'nearest');
colorImReg = imresize(colorIm,scTdf,'nearest');

% topographical distance function, or combination of these
c = tdf(minimaReg,imReg,colorImReg);
for i = 1 : 4
   cFull(:,:,i) = imresize(c(:,:,i),[M N],'bicubic');
end;
c = cFull;

% save the topographical distance function
ver = version;

if strcmp(ver(1:3),'6.5')
    save(name,'c')        
else
    save(name,'c','-V6')    
end;

%----------------------------------------------------------



function [colorIm] = findCol(boundIm)


% how many objercts?
L = max(boundIm(:));
if L <= 4
   colorIm = boundIm; 
    return;
end

try
    L = max(boundIm(:));
    colorIm = zeros(size(boundIm));
    load ball1;se = getball(ball,1,1);
    for i = 1 : L
        regHere = eq(boundIm,i);
        ind = find(regHere);

        % find neighbor values
        neighVal = unique(colorIm(imdilate(regHere,se)));
        % remove the 0 value since it is not a valid value, only means an empty
        % region
        neighVal(neighVal == 0) = [];

        % check the four colors, which of them are NOT in the neighbours
        match = zeros(4,1);        
%        show(regHere,3);colorbar
%        show(colorIm,4);colorbar
        for j = 1 : 4
            indVal = find(neighVal == j);
            if  ~isempty(indVal)
                match(j) = j;
            end;        
        end;
        val = find(match == 0);
        val = val(1);
    %     val
    %     match
    %     show(regHere,2)
    %     pause
    %     
    %     % the values that we have before, to choose four if possible to do
    %     valBefore = unique(colorIm(find(colorIm)));
    %     val = find(match == 0);
    %     for j = 1 : size(val,1)
    %         valHere = val(j)
    %         if isempty(find(valBefore == valHere))
    %             val = [];
    %             val = valHere;
    %             break;
    %         end;
    %           
    %     end;
        colorIm(ind) = val;
    %    show(colorIm,3);colorbar
    %    pause
    end;
catch
    colorIm = zeros(size(boundIm));
    for i = 1 : max(boundIm(:))
        i
        regHere = eq(boundIm,i);
        show(regHere,2);colorbar
        show(colorIm,3);colorbar        
        v = input('Choose color 1 --> 4');        
        colorIm(regHere) = v; 
    end;
end;
% show(boundIm,2);colorbar
%show(colorIm,3);colorbar
%pause

%---------------------------------------------------


function [gamma,epsPCLS,lambda,vis,sigma] = getprm(solver,prm)


    
try
    gamma = prm.gamma;
catch
    gamma = 0.05;
end;
    
try
    epsPCLS = prm.epsPCLS;
catch
    epsPCLS = 2*gamma;
end;
     
try
    vis = prm.vis;
catch
    vis =inf;
end;
        

try
    lambda = prm.lambda;
catch
    lambda = 0.1;
end;

try
    sigma = prm.sigma;
catch
    sigma = 0.01;
end;

