%Chapter 2, Fig. 2.13 and Fig. 2.14
%Eric Dubois, updated 2018-11-08
%Perspective plot of Laplacian of Gaussian function: impulse response and
%frequency response

clear all, close all
r0 = .0025;
%Frequency response
u = -256:16:256; v = u;
[U,V] = meshgrid(u,v);
H = (U.^2 + V.^2).*exp(-2*(pi^2)*(r0^2)*(U.^2 + V.^2));
c = 1/max(max(H));
H = c*H;
figure

mesh(U,V,H)            % generate the perspective plot
colormap([0 0 0]);     % use black only
set(gca,'ydir','reverse'); 
xlabel('\itu \rm(c/ph)'), ylabel('\itv \rm(c/ph)');
axis([-300 300 -300 300 0 1]);

set(gca,'FontName','times')          %set the plot font to Times
set(gca,'FontSize',8)                %set the plot size to 8pt
set(gcf,'units','centimeters','position',[0.5,0.5,15,15]);
set(gcf,'Color',[1 1 1]);

%impulse response
x=-.01:.00075:.01; y=x;        %x and y values between -1 and +1 spaced by 0.025
[X,Y]=meshgrid(x,y);   % generate a 2D grid of xy values
h = -(X.^2 + Y.^2 -2*r0^2).*exp(-(X.^2 + Y.^2)/(2*r0^2));   % generate the Laplacian of Gaussian function on the grid
h = h*c/(2*pi*r0^2)^3;
figure

mesh(X,Y,h)            % generate the perspective plot
colormap([0 0 0]);     % use black only
xlabel('\itx \rm(ph)'), ylabel('\ity \rm(ph)');
set(gca,'ydir','reverse'); 
axis([-.01 .01 -.01 .01 -1E4 8E4]);

set(gca,'FontName','times')          %set the plot font to Times
set(gca,'FontSize',8)                %set the plot size to 8pt
set(gcf,'units','centimeters','position',[0.5,0.5,15,15]);
set(gcf,'Color',[1 1 1]);

% Filter grayscale Barbara image with such a filter
%Create impulse response on a square lattice with spacing 1/512 ph. Scale
%the above response by 1/(512)^2. Make a 21 x 21 filter
xd = -10/512:1/512:10/512; yd = xd;
[Xd,Yd] = meshgrid(xd,yd);
hd = (Xd.^2 + Yd.^2 - 2*r0^2).*exp(-(Xd.^2 + Yd.^2)/(2*r0^2)); 
hd = hd*c/(2*pi*r0^2)^3/(512^2);
IMG = im2double(imread('barb512_gray.tif'));
IMG_Filt = imfilter(IMG,hd,'symmetric','same');
figure,imshow(IMG_Filt+.5)