level set algorithm

close all;

Img=minion ;
ask="Do you want to add noise ?(y/n) : ";
noise_check=input(ask,"s");
if noise_check=='y'
    noise_weight="Enter the noise weight : ";
    noise_w=input(noise_weight);
    Img= imnoise(Img,'gaussian',noise_w);
end
%{
ask="Enter a letter : ";
letter=input(ask,"s");
%}
Img=double(Img(:,:,1));
original=Img;

%% parameter setting
timestep=5;  % time step
mu=0.2/timestep;  % coefficient of the distance regularization term R(phi)
iter_inner=2;
iter_outer=50;
lambda=5; % coefficient of the weighted length term L(phi)
alfa=1.5;  % coefficient of the weighted area term A(phi)
epsilon=1.5; % papramater that specifies the width of the DiracDelta function
sigma=1.5;     % scale parameter in Gaussian kernel
G=fspecial('gaussian',15,sigma);
Img_smooth=conv2(Img,G,'same');  % smooth image by Gaussiin convolution
[Ix,Iy]=gradient(Img_smooth);
f=Ix.^2+Iy.^2;
g=1./(1+f);  % edge indicator function.

% initialize LSF as binary step function
c0=2;
initialLSF=c0*ones(size(Img));
figure
imshow(minion);
[x,y]=ginput(2);
close
x=round(x);
y=round(y);
% generate the initial region R0 as a rectangle
initialLSF(y(1):y(2), x(1):x(2))=-c0;  
phi=initialLSF;

figure(1);
mesh(-phi);   % for a better view, the LSF is displayed upside down
hold on;  contour(phi, [0,0], 'r','LineWidth',2);
title('Initial level set function');
view([-80 35]);

figure(2);
imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on;  contour(phi, [0,0], 'r');
title('Initial zero level contour');
pause(0.5);

potential=2;  
if potential ==1
    potentialFunction = 'single-well';  % use single well potential p1(s)=0.5*(s-1)^2, which is good for region-based model 
elseif potential == 2
    potentialFunction = 'double-well';  % use double-well potential in Eq. (16), which is good for both edge and region based models
else
    potentialFunction = 'double-well';  % default choice of potential function
end


% start level set evolution
for n=1:iter_outer
    phi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction);
    if mod(n,2)==0
        figure(2);
        imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on;  contour(phi, [0,0], 'r');
    end
end

% refine the zero level contour by further level set evolution with alfa=0
alfa=0;
iter_refine = 10;
phi = drlse_edge(phi, g, lambda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction);
%{
if letter=='a'
    finalLSF_high=phi;
    finalLSF= finalLSF_high;
else
    finalLSF=phi;
end
%}
finalLSF=phi;
figure(2);
imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on;  contour(phi, [0,0], 'r');
hold on;  contour(phi, [0,0], 'r');
str=['Final zero level contour, 150 iterations'];
title(str);

pause(1);
figure;
mesh(-finalLSF); % for a better view, the LSF is displayed upside down
hold on;  contour(phi, [0,0], 'r','LineWidth',2);
str=['Final level set function, 150 iterations'];
title(str);
axis on;

if noise_check=='y'
    [ps ,sn]=psnr(finalLSF,finalLSF_no_noise);
    ms=mse(finalLSF_noise,finalLSF_no_noise);
    ps=10*(log((332*208)/ms)/log(10));
    ssi=ssim(finalLSF_noise,finalLSF_no_noise);
    fprintf("PSNR(noise = %f ) : %f\n",noise_w,ps)
    fprintf("SNR(noise = %f ) : %f\n",noise_w,sn)
    fprintf("MSE(noise = %f ) : %f\n",noise_w,ms) %mean square error
    fprintf("SSIM(noise = %f ) : %f\n",noise_w,ssi)%structural similarity index measure 
    end

%{
finalLSF= imresize(finalLSF,size(finalLSF_high));
if letter~='a'
    [ps ,sn]=psnr(finalLSF,finalLSF_high);
    ms=mse(finalLSF,finalLSF_high);
    ps=10*(log((332*208)/ms)/log(10));
    ssi=ssim(finalLSF,finalLSF_high);
    fprintf("PSNR(low_contrast) : %f\n",ps)
    fprintf("SNR(low_contrast) : %f\n",sn)
    fprintf("MSE(low_contrast) : %f\n",ms) %mean square error
    fprintf("SSIM(low_contrast) : %f\n",ssi)%structural similarity index measure 
end
%}