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

%Initialization
X=100;
Y=100;
radius=25;
center=[X/2,Y/2];
scaleIso=95;
stopIso=1e-5;
scaleNC=35;
stopNC=1e-2;
tiles=50;
shotMult=5000;
minNoise=0;
maxNoise=.2;

%Generate image
img=zeros(X,Y);
[y x]=meshgrid(1:X,1:Y);
index=find(sqrt((x-center(1)).^2 + (y-center(2)).^2)<radius);
img(index)=1;

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

    %Additive Gaussian noise    
    noiseImg=noiseVar*randn(X,Y)+img;
    [imgmasks,imgMarkup,segOutline]=imgsegment(noiseImg,scaleIso,stopIso);
    countAddIso(count)=max(max(imgmasks));    
    [imgmasks, imgMarkup, segOutline]= ...
        imgsegment(noiseImg,scaleNC,stopNC,1,0,1);
    countAddNC(count)=max(max(imgmasks));    
    
    %Multiplicative Gaussian noise
    noiseImg=img.*(noiseVar*randn(X,Y)+1);
    [imgmasks,imgMarkup,segOutline]=imgsegment(noiseImg,scaleIso,stopIso);
    countMultIso(count)=max(max(imgmasks));        
    [imgmasks, imgMarkup, segOutline]= ...
        imgsegment(noiseImg, scaleNC, stopNC, 1, 0,1);
    countMultNC(count)=max(max(imgmasks));    
    
    %Shot noise
    numbShots=round(shotMult*noiseVar);
    tmp=randperm(X*Y);
    img2=img;
    img2(tmp(1:numbShots))=1;
    [imgmasks,imgMarkup,segOutline]=imgsegment(img2,scaleIso,stopIso);
    countShotIso(count)=max(max(imgmasks));            
    [imgmasks,imgMarkup,segOutline]= ...
        imgsegment(img2,scaleNC,stopNC,1,0,1);
    countShotNC(count)=max(max(imgmasks));                
    [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
subplot(3,1,1)
plot(linspace(minNoise,maxNoise,tiles),countAddIso(1:tiles),'k')
hold on
plot(linspace(minNoise,maxNoise,tiles),countAddNC(1:tiles),'k--')
xlabel('Noise variance')
ylabel('Number of segments found')
title('Additive noise')
axis tight
axis([0 noiseVar 0 ...
        max(max(countAddIso(1:tiles)),max(countAddNC(1:tiles)))])

subplot(3,1,2)
plot(linspace(minNoise,maxNoise,tiles),countMultIso(1:tiles),'k')
hold on
plot(linspace(minNoise,maxNoise,tiles),countMultNC(1:tiles),'k--')
xlabel('Noise variance')
ylabel('Number of segments found')
title('Multiplicative noise')
axis tight
axis([0 noiseVar 0 ...
        max(max(countMultIso(1:tiles)),max(countMultNC(1:tiles)))])

subplot(3,1,3)
plot(round(shotMult*(linspace(minNoise,maxNoise,tiles))),countShotIso(1:tiles),'k')
hold on
plot(round(shotMult*(linspace(minNoise,maxNoise,tiles))),countShotNC(1:tiles),'k--')
xlabel('Number of shots')
ylabel('Number of segments found')
title('Shot noise')
axis tight
axis([0 round(shotMult*noiseVar) 0 ...
        max(max(countShotIso(1:tiles)),max(countShotNC(1:tiles)))])