%Script used to generate segmentations using different methods for a 
%simple image.
%
%
%6/5/03 - Leo Grady

% Copyright (C) 2002, 2003 Leo Grady <lgrady@cns.bu.edu>
%   Computer Vision and Computational Neuroscience Lab
%   Department of Cognitive and Neural Systems
%   Boston University
%   Boston, MA  02215
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
%
% Date - $Id: segmentationCompareDemo.m,v 1.4 2003/08/21 17:29:30 lgrady Exp $
%========================================================================%
%global LPIRVELI LMEORE
clear all
close all

addpath ../../graphAnalysisToolbox-1.0/

global LPIRVELI LMEORE

i=0;

%global LPIRVELI LMEORE

%Initialization
stopIso=1;
stopNcuts=5e-2;
scaleIso=100;
scaleNcuts=100;

%Read image
%img=im2double(imread('./../images/blood1.tif'));
%image=im2double(imread('/home/guests/ajujunas/SpectralProject/lenna.gif'));
%image=im2double(imread('/home/guests/ajujunas/SpectralProject/Forest 118_1830.JPG'));
image=im2double(imread('/home/guests/ajujunas/SpectralProject/phot-02a-98-hires.jpg'));
%Resize to make script run faster
[X Y Z]=size(image);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%How many times we reduce the picture at the beginning?
nBegin=64;
%How many pictures do we need?
nPictures=3;
%What is the last resizer?
nLast=128;
%The rseize number at each step:
resizemult=1;
if (nPictures>1)
    resizemult=(nBegin/nLast)^(1/(nPictures-1));
end %end if
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
image1=imresize(image,[floor(X/nBegin),floor(Y/nBegin)]);
[X Y Z]=size(image1);
%From what level begin the partitioning?
nLevel=0; % nLevel must less than nPictures
for ii=nLevel:(nPictures-1)
    img=imresize(image1,[floor(X*(resizemult^ii)),floor(Y*(resizemult^ii))]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%img=imresize(img,[floor(X/nBegin),floor(Y/nBegin)]);


%Perform isoperimetric segmentation
% tic
% [imgmasksIso,imgMarkupIso,segOutlineIso]= ...
%     imgsegment(img,scaleIso,stopIso);
% isoTime=toc

%Perform spectral segmentation
tic
[imgmasksSpec,imgMarkupSpec,segOutlineSpec]= ...
    imgsegment(img,scaleIso,stopIso,1,1,1);
specTime=toc

%Perform Ncuts segmentation
tic
[imgmasksNC,imgMarkupNC,segOutlineNC]= ...
    imgsegment(img,scaleNcuts,stopNcuts,1,0,1);
ncutsTime=toc


% %Display
% figure
% subplot(2,2,1)
% imagesc(img)
% axis equal
% axis tight
% axis off
% colormap(gray)
% title('Original image')
% 
% subplot(2,2,2)
% imagesc(imgMarkupIso)
% axis equal
% axis tight
% axis off
% colormap(gray)
% title('Isoperimetric segmentation')
% 
% subplot(2,2,3)
% imagesc(imgMarkupSpec)
% axis equal
% axis tight
% axis off
% colormap(gray)
% title('Spectral segmentation')
% 
% subplot(2,2,4)
% imagesc(imgMarkupNC)
% axis equal
% axis tight
% axis off
% colormap(gray)
% title('NCuts segmentation')
i=i+1;
options.tol=10^-2;
%tic
t=cputime;
[v,dL,flag]=eigs(LPIRVELI,1,'la',options);
%%%%%%%%%%%%%%%%%%%
%d=eigs(L,1);
%%%%%%%%%%%%%%%%%%%
%timeLLarg(i)=toc;
timeLLarg(i)=cputime-t;
ns=size(LPIRVELI);
sizeL(i)=ns(1);
LambLargL(i)=dL;

LL2=(LMEORE+LMEORE')/2;

%tic
t=cputime;
[v,dL2,flag]=eigs(LL2,1,'la',options);
%%%%%%%%%%%%%%%%%%%
%d=eigs(LL2,1);
%%%%%%%%%%%%%%%%%%%
%timeL2Larg(i)=toc;
timeL2Larg(i)=cputime-t;
LambLargL2(i)=dL2;


%tic
t=cputime;
dLSm=eigs(LPIRVELI,2,10^-7);
%timeLSmall2(i)=toc;
timeLSmall2(i)=cputime-t;
LambSmalL(i)=dLSm(1);


% Calculating of a few first clustered Eigenvalues

clusterEigsL=eigs(LPIRVELI,11,10^-7);
lambSmal_2_L(i)=clusterEigsL(10);
lambSmal_3_L(i)=clusterEigsL(9);
lambSmal_4_L(i)=clusterEigsL(8);
lambSmal_5_L(i)=clusterEigsL(7);
lambSmal_6_L(i)=clusterEigsL(6);
lambSmal_7_L(i)=clusterEigsL(5);
lambSmal_8_L(i)=clusterEigsL(4);
lambSmal_9_L(i)=clusterEigsL(3);
lambSmal_10_L(i)=clusterEigsL(2);
lambSmal_11_L(i)=clusterEigsL(1);


%tic
t=cputime;
dL2Sm=eigs(LL2,2,10^-7);
%timeL2Small2(i)=toc;
timeL2Small2(i)=cputime-t;
LambSmalL2(i)=dL2Sm(1);


% Calculating of a few first clustered Eigenvalues

clusterEigsL2=eigs(LL2,11,10^-7);
lambSmal_2_L2(i)=clusterEigsL2(10);
lambSmal_3_L2(i)=clusterEigsL2(9);
lambSmal_4_L2(i)=clusterEigsL2(8);
lambSmal_5_L2(i)=clusterEigsL2(7);
lambSmal_6_L2(i)=clusterEigsL2(6);
lambSmal_7_L2(i)=clusterEigsL2(5);
lambSmal_8_L2(i)=clusterEigsL2(4);
lambSmal_9_L2(i)=clusterEigsL2(3);
lambSmal_10_L2(i)=clusterEigsL2(2);
lambSmal_11_L2(i)=clusterEigsL2(1);


CondL(i)=LambLargL(i)/LambSmalL(i);
CondL2(i)=LambLargL2(i)/LambSmalL2(i);

%%% From here we'll make comments
%Display

figure
subplot(1,5,1)
imagesc(img)
axis equal
axis tight
axis off
colormap(gray)
%title('Orig.image')
%title('Size of Image - ''Pixels')

subplot(1,5,2)
imagesc(imgMarkupSpec)
axis equal
axis tight
axis off
colormap(gray)
%title('Spec.segment')

subplot(1,5,3)
imagesc(imgMarkupNC)
axis equal
axis tight
axis off
colormap(gray)
%title('NCut.segment')

subplot(1,5,4)
imagesc(segOutlineSpec)
axis equal
axis tight
axis off
colormap(gray)
%title('Spect.line')

subplot(1,5,5)
imagesc(segOutlineNC)
axis equal
axis tight
axis off
colormap(gray)
%title('NCuts.line')
%hh=figure;
%print -depsc2 hh; 
%%%% Before here comments to all lines

end %End of cycle


%CondLlog=log(CondL);
%CondL2log=log(CondL2);
%sizeLlog=log(sizeL);

%%%% Plotting of Cond. Numbers

figure
%subplot(2,1,1)
title('Condition Numbers of Matrices L and L2 in Subspaces (lin. scal)')
plot(sizeL,CondL,'r',sizeL,CondL2,'--g')
% plot(sizeL,CondL);
% 
% figure
% plot(sizeL,CondL2);

%subplot(2,1,2)
%title('Condition Numbers of Matrices L and L2 in Subspaces (log. scal)')
%plot(sizeLlog,CondLlog,'r*',sizeLlog,CondL2log,'g-.')
%plot(sizeL,CondL,'r',sizeL,CondL2,'--g')

%%%% Plotting of Elapsed Times

%figure
%subplot(2,1,1)
%title('Time Lambda Larg')
%plot(sizeL,timeLLarg,'r*',sizeL,timeL2Larg,'g-.')
% plot(sizeL,CondL);
% 
% figure
% plot(sizeL,CondL2);

figure
%subplot(2,1,2)
title('Time Lambda Small')
%plot(sizeLlog,CondLlog,'r*',sizeLlog,CondL2log,'g-.')
plot(sizeL,timeLSmall2,'r',sizeL,timeL2Small2,'--g')

%%%% Plotting of the small eigenvalues

figure
subplot(2,1,1)
title('A few smal eigenvalues of matrix L')
plot(sizeL,lambSmal_2_L,'y',sizeL,lambSmal_3_L,'m',sizeL,lambSmal_4_L,'c',sizeL,lambSmal_5_L,'r',sizeL,lambSmal_6_L,'g',sizeL,lambSmal_7_L,'b',sizeL,lambSmal_8_L,'k',sizeL,lambSmal_9_L,'y*',sizeL,lambSmal_10_L,'c+',sizeL,lambSmal_11_L,'mo');


subplot(2,1,2)
title('A few smal eigenvalues of matrix L2')
plot(sizeL,lambSmal_2_L2,'y',sizeL,lambSmal_3_L2,'m',sizeL,lambSmal_4_L2,'c',sizeL,lambSmal_5_L2,'r',sizeL,lambSmal_6_L2,'g',sizeL,lambSmal_7_L2,'b',sizeL,lambSmal_8_L2,'k',sizeL,lambSmal_9_L2,'y*',sizeL,lambSmal_10_L2,'c+',sizeL,lambSmal_11_L2,'mo');


% 
% figure
% plot(sizeL,timeLLarg);
% 
% figure
% plot(sizeL,timeL2Larg);
% 
% figure
% plot(sizeL,timeLSmall2);
% 
% figure
% plot(sizeL,timeL2Small2);