imbase = ['../text/images/'];
% Get some images
bi = imread([imbase,'circles.tif']);     % Binary
gs = imread([imbase,'cameraman.tif']);   % Grayscale
co = imread([imbase,'flowers.tif']);     % RGB color
[in inmap] = imread([imbase,'emu.tif']); % Indexed color

% 1D Gausian
n = (-50:50)';
sigma1 = 5;
sigma2 = 10;
h1 = exp(-n.^2/(2*sigma1^2));h1 = h1/sum(h1);
h2 = exp(-n.^2/(2*sigma2^2));h2 = h2/sum(h2);
figure(1);clf
plot(n,h1,'-',n,h2,'--');
legend('sigma = 5','sigma = 10');
print -djpeg gaussian1d.jpg
[sum(h1) sum(h2)]

% Use of meshgrid
n = -1:1;
m = -1:1;
[M N] = meshgrid(m,n)

% 2D Gausian (simple case);
n = (-40:40)';
m = n;
[M N] = meshgrid(m,n);

h1 = exp(-1/(2*sigma1^2)*(M.^2 + N.^2));
c = sum(sum(h1)); h1 = h1/c;
h2 = exp(-1/(2*sigma2^2)*(M.^2 + N.^2));
c = sum(sum(h2)); h2 = h2/c;

figure(1);clf;colormap jet
subplot(221);mesh(m,n,h1)
subplot(222);mesh(m,n,h2)
subplot(223);imagesc(m,n,h1);axis image
subplot(224);imagesc(m,n,h2);axis image
print -djpeg gaussian2d.jpg

% Box filters
h1 = fspecial('average',[3 3]);
h2 = fspecial('average',[5 5]);
h3 = fspecial('average',[13 1]);
h4 = fspecial('average',[1 17]);

% Cameraman
f = double(gs);
% Dummy image
f = [10*ones(128) 100*ones(128); ...
  50*ones(128) 200*ones(128)];

f1 = filter2(h1,f,'same');
f2 = filter2(h2,f,'same');
f3 = filter2(h3,f,'same');
f4 = filter2(h4,f,'same');

figure(1);clf;colormap gray
imagesc(f);axis image

figure(2);clf;colormap(gray)
subplot(221);imagesc(f1);axis square
subplot(222);imagesc(f2);axis square
subplot(223);imagesc(f3);axis square
subplot(224);imagesc(f4);axis square

% Gaussian filters
h1 = fspecial('gaussian',20,1);
h2 = fspecial('gaussian',20,3);
h3 = fspecial('gaussian',20,7);
h4 = fspecial('gaussian',20,10);

figure(1);clf;colormap hsv
subplot(221);imagesc(h1);axis image
title('\sigma = 1');
subplot(222);imagesc(h2);axis image
title('\sigma = 3');
subplot(223);imagesc(h3);axis image
title('\sigma = 7');
subplot(224);imagesc(h4);axis image
title('\sigma = 10');

% Cameraman
f = double(gs);
% Dummy image
f = [10*ones(128) 100*ones(128); ...
  50*ones(128) 200*ones(128)];

f1 = filter2(h1,f,'same');
f2 = filter2(h2,f,'same');
f3 = filter2(h3,f,'same');
f4 = filter2(h4,f,'same');

figure(1);clf;colormap gray
imagesc(f);axis image

figure(2);clf;colormap(gray)
subplot(221);imagesc(f1);axis square
subplot(222);imagesc(f2);axis square
subplot(223);imagesc(f3);axis square
subplot(224);imagesc(f4);axis square

% Something cool that works only with Gaussians
f = double(gs);
f = [10*ones(128) 100*ones(128); ...
  50*ones(128) 200*ones(128)];
h1 = fspecial('gaussian',41,3);
h2 = fspecial('gaussian',41,3.1);
f1 = filter2(h1,f,'same');
f2 = filter2(h2,f,'same');
figure(1);clf;colormap gray
imagesc(f);axis image
figure(2);clf;colormap(gray)
imagesc(f1-f2);axis square;colorbar