%Script used to perform noise analysis and comparison between the 
%isoperimetric and Ncuts algorithms
%
%Note: Code involves randomness, so results may be slightly different from
%       those in thesis.
%
%
%6/6/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: fig4.m,v 1.1 2006/02/19 21:06:16 jonnyreb Exp $
%========================================================================%

clear
close all

%Initial parameters
tiles=5;
minNoise=.01;
maxNoise=.2;
scale=95; %Different than standard to accomodate Ncuts
stop=1e-5;
ncutsScale=35;
radius=0;
ncutsStop=5e-2; %Different since ncuts criterion is different than iso

%Read image
img=im2double(imread('blood1.tif'));
[X Y Z]=size(img);
img=imresize(img,[floor(X/2),floor(Y/2)]);
[X Y Z]=size(img);

%Initialize tiled images
tileFig=zeros((tiles+1)*X,9*Y);

%Initial segmentation
[imgmasks,imgMarkup,segOutline]=imgsegment(img,scale,stop);
[segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
        segOutline(X,1:Y)]=deal(0);
tileFig(1:X,1:Y)=img;
tileFig(1:X,(Y+1):(2*Y))=segOutline;
[imgmasks,imgMarkup,segOutline]= ...
    imgsegment(img,ncutsScale,ncutsStop,1,0,1);
[segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
        segOutline(X,1:Y)]=deal(0);
tileFig(1:X,(2*Y+1):(3*Y))=segOutline;
tileFig(1:X,(3*Y+1):9*Y)=[tileFig(1:X,1:3*Y),tileFig(1:X,1:3*Y)];

%Progressively add noise
count=0;
for noiseVar=linspace(minNoise,maxNoise,tiles)
    count=count+1;

    %Additive Gaussian noise    
    noiseImg=noiseVar*randn(X,Y)+img;
    noiseImg2=noiseImg-min(min(noiseImg));
    noiseImg2=noiseImg2./max(max(noiseImg2));
    tileFig((count*X+1):(count+1)*X,1:Y)=noiseImg2;
    [imgmasks,imgMarkup,segOutline]=imgsegment(noiseImg,scale,stop);
    [segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
            segOutline(X,1:Y)]=deal(0);
    tileFig((count*X+1):(count+1)*X,(Y+1):2*Y)=segOutline;
    [imgmasks, imgMarkup, segOutline]= ...
        imgsegment(noiseImg,ncutsScale,ncutsStop,1,0,1);
    [segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
            segOutline(X,1:Y)]=deal(0);
    tileFig((count*X+1):(count+1)*X,(2*Y+1):3*Y)=segOutline;

    %Multiplicative Gaussian noise
    noiseImg=img.*(noiseVar*randn(X,Y)+1);
    noiseImg2=noiseImg-min(min(noiseImg));
    noiseImg2=noiseImg2./max(max(noiseImg2));    
    tileFig((count*X+1):(count+1)*X,(3*Y+1):4*Y)=noiseImg2;
    [imgmasks,imgMarkup,segOutline]=imgsegment(noiseImg,scale,stop);
    [segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
            segOutline(X,1:Y)]=deal(0);
    tileFig((count*X+1):(count+1)*X,(4*Y+1):5*Y)=segOutline;
    [imgmasks, imgMarkup, segOutline]= ...
        imgsegment(noiseImg,ncutsScale,ncutsStop,1,0,1);
    [segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
            segOutline(X,1:Y)]=deal(0);
    tileFig((count*X+1):(count+1)*X,(5*Y+1):6*Y)=segOutline;

    %Shot noise
    numbShots=round(1000*noiseVar);
    tmp=randperm(X*Y);
    noiseImg=img;
    noiseImg(tmp(1:numbShots))=1;
    noiseImg2=noiseImg-min(min(noiseImg));
    noiseImg2=noiseImg2./max(max(noiseImg2));    
    tileFig((count*X+1):(count+1)*X,(6*Y+1):7*Y)=noiseImg2;    
    [imgmasks,imgMarkup,segOutline]=imgsegment(noiseImg,scale,stop);
    [segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
            segOutline(X,1:Y)]=deal(0);
    tileFig((count*X+1):(count+1)*X,(7*Y+1):8*Y)=segOutline;
    [imgmasks,imgMarkup,segOutline]= ...
        imgsegment(noiseImg,ncutsScale,ncutsStop,1,0,1);
    [segOutline(1:X,1),segOutline(1,1:Y),segOutline(1:X,Y), ...
            segOutline(X,1:Y)]=deal(0);
    tileFig((count*X+1):(count+1)*X,(8*Y+1):9*Y)=segOutline;
end

%Display
figure
imagesc(tileFig)
colormap(gray)
axis equal
axis tight
axis off