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

% Why we need smoothing for edge detection

f = zeros(128);
f(:,64:end) = 1;
f = f+.5*randn(128);
h1 = [-1 0 1];
h2 = [-1 0 1 ; -1 0 1; -1 0 1];
h2 = [-1 0 1 ; -1 0 1; -1 0 1];
h3 = [-1 0 1 ; -2 0 2; -1 0 1];
df1 = filter2(h1,f,'same');
df2 = filter2(h2,f,'same');
df3 = filter2(h3,f,'same');
figure(1);clf;colormap gray
subplot(221);imagesc(f);axis image
title('Noisy image with edge')
subplot(222);imagesc(abs(df1));axis image
title('No smoothing')
subplot(223);imagesc(abs(df2));axis image
title('Prewitt')
subplot(224);imagesc(abs(df3));axis image
title('Sobel')

% Laplacian of a Gaussian
x = linspace(-10,10,101)';
y = x;
[X Y] = meshgrid(x,y);
sig = 2;
tmp1 = (X.^2/(sig^4) + Y.^2/(sig^4) - 2/(sig^2));
tmp2 = exp(-(X.^2+Y.^2)/(2*sig*sig));
LapGau = tmp1.*tmp2;
figure(1);clf;colormap hsv
surf(x,y,-LapGau);shading interp
print -djpeg LapGau_surf.jpg

figure(1);clf
plot(x,-LapGau(:,50))
title('Slice of Negative of Laplacian of Gaussian')
print -djpeg LapGau_slice.jpg

% Some edge filters
% Dummy image
f = [10*ones(128) 100*ones(128); ...
  50*ones(128) 200*ones(128)];
fn = f+20*randn(size(f));

% Sobel filters
s_h  = fspecial('sobel'); s_v = s_h';
dfh = filter2(s_h,f);
dfv = filter2(s_v,f);
df  = sqrt(dfh.^2 + dfv.^2);

dfnh = filter2(s_h,fn);
dfnv = filter2(s_v,fn);
dfn  = sqrt(dfnh.^2 + dfnv.^2);

% Laplacian
lap = 0.25*[1 2 1 ; 2 -12 2 ;1 2 1 ];
dflap = filter2(lap,f);
dfnlap = filter2(lap,fn);

% Laplacian of Gaussian
log = fspecial('log',25,3);
dflog = filter2(log,f);
dfnlog = filter2(log,fn);