%% compare Gaussian and Poisson noise

% im = double(phantom(256));
% im = im - 2*(min(min(im)));
N = 256;
im = phantom(N);
im = im - min(im(:));

figure(5);
subplot(2,2,1);
imagesc(im);colorbar;title('original f');
subplot(2,2,2);
scale = 1e12;
yn12  = scale*imnoise(im/scale,'poisson');
imagesc(yn12);colorbar;title('Poisson noise; \sigma = f');
subplot(2,2,3);
scale = 1e10;
yn10  = scale*imnoise(im/scale,'poisson');
imagesc(yn10);colorbar;title('Poisson noise; \sigma = 10% of f');
subplot(2,2,4);
scale = 1e8;
yn8  = scale*imnoise(im/scale,'poisson');
imagesc(yn8);colorbar;title('Poisson noise; \sigma =  1% of f');
%subplot(2,2,3);
%imagesc(yn-im);colorbar;title('Poisson noise');

%% compare to Gaussian noise
figure(6);

subplot(2,2,1);
imagesc(im);colorbar;title('original f');
subplot(2,2,2);
yg12 = im + sqrt(im).*randn(size(im));
imagesc(yg12);colorbar;title('Gaussian noise; \sigma= 100% of f');
subplot(2,2,3);
yg10 = im + 0.1*sqrt(im).*randn(size(im));
imagesc(yg10);colorbar;title('Gaussian noise; \sigma= 10% of f ');
subplot(2,2,4);
yg8 = im + 0.01*sqrt(im).*randn(size(im));
imagesc(yg8);colorbar;title('Gaussian noise; \sigma = 1% of f');

%% compare noise
ng8 = yg8 - im;
np8 = yn8 - im;

figure(7);clf;
subplot(2,2,1); imagesc(ng8);colorbar; title('Additive Gaussian noise 1%');
subplot(2,2,2); imagesc(np8);colorbar;  title('Additive Poisson noise 1%');


ng12 = yg12 - im;
np12 = yn12 - im;

subplot(2,2,3); imagesc(ng12);colorbar;  title('Additive Gaussian noise 100%');
subplot(2,2,4); imagesc(np12);colorbar;  title('Additive Poisson noise 100%');

%% Historgrams?
%Histograms are not too helpful, because variance is spatially varying...
figure(8);
subplot(2,2,1); histogram(ng8); title('Additive Gaussian noise 1%');
subplot(2,2,2); histogram(np8); title('Additive Poisson noise 1%');
subplot(2,2,3); histogram(ng12); title('Additive Gaussian noise 100%');
subplot(2,2,4); histogram(np12); title('Additive Poisson noise 100%');