clc;clear all;close all;
%------------------------------------
% This Matlab program demomstrates the level set algorithm in paper:
%   Kaihua Zhang, Lei Zhang,and Su Zhang "A VARIATIONAL MULTIPHASE LEVEL SET APPROACH TO SIMULTANEOUS SEGMENTATION AND BIAS CORRECTION" 
%    ICIP, 2010
% Author: Kaihua Zhang
% E-mail: cskhzhang@comp.polyu.edu.hk, cslzhang@comp.polyu.edu.hk  
% http://www4.comp.polyu.edu.hk/~cskhzhang/
% http://www4.comp.polyu.edu.hk/~cslzhang/
%  Notes:
%   1. Some parameters may need to be modified for different types of images. Please contact the author if any problem regarding the choice of parameters.
%   2. Intial contour should be set properly.
% Date: 19/6/2010
%------------------------------------Set path
fp = mfilename('fullpath');
rootdir = fileparts(fp);
p = fullfile(rootdir,'image');
res = fullfile(rootdir,'phi_initialization');
addpath(rootdir,p,res);
%------------------------------------
I = imread('3Tbrain2.png'); % for 3t MRI
%   I = imread('brain9.jpg');%for 7t MRI
I = double(I(:,:,1));
[row,col] = size(I);

sigma = 4.5;
K_R = 1/9*ones(3);% The constant kernel for regularizing the level set function after each iteration
K = ones(4*sigma+1);%The constant kerenl for SVMLS
%initialization
dim = 4; % the number of segmentation regions

phi_flag = 1;%flag of different intializations: 1;3t MRI,aross the image, 2: outside the image, 3: 7t MRI.Different initial contours.
switch phi_flag
    case 1;
        load phi_cross_3t.mat
    case 2;
        load phi_out_3t.mat;
    case 3; 
        load phi_cross_7t.mat;
end
phi1= phi_1;
phi2= phi_2;
% Initialize b,s,c.   
b(1:row,1:col) = abs(random('norm',0,1,row,col));
for i = 1:dim
s(1:row,1:col,i) = abs(random('norm',0,1,row,col));
end
s(:,:,:) = 1/sqrt(2*pi);
b(:,:) = 1;
for i = 1:dim
c(1:row,1:col,i) = i*(max(I(:))-min(I(:)))./dim + min(I(:));
end
%---------------------------------
Iternum =15; %Iterations
epsilon = 1;
timestep = .45;

for i = 1:Iternum
    u = compute_u(phi1,epsilon,phi2);% the membership function, i.e., u(:,:,1)=M1,...u(:,:,4)=M4 as Eq.(11)
    c = compute_c(I,K,u,b,s);% ci in Eq.(12)
    b = compute_b(I,K,u,c,s);% b in Eq.(12)
    s = compute_s(I,b,K,c,u);% the variance of corresponding region. see the sigma in Eq.(12)
    d = computer_d(I,K,s,b,c);%the d in Eq.(10)
    [phi1,phi2] = evolution(phi1,d,epsilon,timestep,phi2);%Level set evolution equation
    %-------------Regularization procedure
    phi1 = conv2(phi1,K_R,'same');
    phi2 = conv2(phi2,K_R,'same');
    %------------------------------------
   
    pause(0.0001);
    imagesc(I,[0 255]);colormap(gray)
    hold on;
    contour(phi1,[0 0],'r');
    if dim == 4
       contour(phi2,[0 0],'b');
    end   
    iterNum=[num2str(i), ' iterations'];  
    title(iterNum);hold off;        
end

close all;
%------------------------J shows the final segmentation results by regions
J = zeros(row,col);
for i = 1:dim
    J = J + u(:,:,i).*c(:,:,i);
end
%------------------------
if phi_flag == 3;
    M = ones(row,col);
else
   load M;% mask for eliminating the trivial contours outside the object
end
%-------------------------
figure;subplot(1,4,1);imagesc(I,[0 255]);axis off;colormap(gray);
hold on;
contour(phi_1,[0 0],'r');contour(phi_2,[0 0],'b');
hold off;
subplot(1,4,2);imagesc(b.*M);axis off;
subplot(1,4,3);imagesc(J.*M);axis off;
subplot(1,4,4);imagesc(uint8(I.*M./b));axis off;
figure;
subplot(1,2,1);imhist(uint8(I));axis([0 255 0 1000]);title('Hist of Original Image');
subplot(1,2,2);imhist(uint8(I.*M./b));axis([0 255 0 1000]);title('Hist of SVMLS');colormap(gray);

