%Chapter 10, Example 10.2, Fig 10.5-10.7
%Eric Dubois, updated 2018-12-18
%Example for Gaussian bandstop filter design and application to Barbara
close all; clear all;
A = im2double(imread('barb512_gray.tif'));
N = 512;
imshow(A)
L = 12;  %chosen filter half size
x=-L:1:L; y=x;         %x and y values between -L and +L spaced by 1
[X,Y]=meshgrid(x,y);   % generate a 2D grid of xy values
uc = 15/N;  %3dB cutoff in samples/px
r=(.1325/uc);  
h1=exp(-(X.^2+Y.^2)/(2*r^2)); % generate the Gaussian function on the grid
c = 1/(ones(1,2*L+1)*h1*ones(2*L+1,1)); % normalize for unit DC gain
c_approx = 1/(2*pi*r^2);
%u0 = 156;
%v0=-71;
u0=150;   %alternative frequency 1
v0=-86;
%u0=144;  %alternative frequency 2
%v0=-76;
%generate the bandpass filter
h1 = 2*c*h1.*cos(2*pi*(X*u0+Y*v0)/512);
%generate the bandstop filter
h2 = zeros(2*L+1,2*L+1);
h2(L+1,L+1)=1;
h2 = h2 - h1;
[H2,f1,f2] = freqz2(h2);
[F1,F2]=meshgrid(256*f1,256*f2);
figure;
mesh(F1,F2,H2);
xlabel('Horizontal freq. (c/ph)'), ylabel('Vertical freq. (c/ph)')
zlabel('Magnitude')
colormap([0 0 0]);     % use black only
set(gca,'ydir','reverse');  %make v axis point downward
set(gcf,'Color',[1 1 1]);
v=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.1
figure, [C,h]=contour(F1,F2,H2,v); % generate the contour plot, including values to label contours
axis square;  %make the plot square
xlabel('Horizontal frequency (c/ph)'), ylabel('Vertical frequency (c/ph)')
clabel(C,h)             %label the contours
title('Contour plot of Gabor filter, r2=3.75')
set(gca,'ydir','reverse');  %make v axis point downward
colormap([0 0 0]);     % use black only
%Create the filtered Barbara image
B = conv2(A,h1,'same'); % bandpass image -- just for interest
figure,imshow(B+.5)
C = conv2(A,h2,'same');
figure,imshow(C)
%Closeup of the knee 64 x 64
D = A(290:353,170:233);
%zoom it to 256 x 256
D_up = imresize(D,4,'bicubic');
figure, imshow(D_up);
%Closeup of the notched knee
E = C(290:353,170:233);
%zoom it to 256 x 256
E_up = imresize(E,4,'bicubic');
figure, imshow(E_up);
%Closeup of bandpass filtered knee
F = B(290:353,170:233);
%zoom it to 256 x 256
F_up = imresize(F,4,'bicubic');
figure, imshow(F_up+.5);