%Chapter 10, Example 10.1, Fig 10.2-10.4
%Eric Dubois, updated 2018-12-18
% Example of filter design for lowpass filtering
% Design a low-pass filter on a square lattice with a 3 dB bandwidth
% of 1/8X c/ph (1/8 c/px) in both the horizontal and the vertical directions (3 dB bandwidth
% corresponds to an attenuation of 0.707).
% Use the moving average filter, the Gaussian filter, the window method.
% Frequencies will be measured in c/px.
% Apply the filters to the barbara image.
%
%
clear all, close all

% Moving average filter
% To achieve the required attenuation, we require $L=2$, i.e., a $5 \times 5$ 
% moving average filter. A red circle of radius 1/8 c/px is plotted as well. We
% see that the attenuation of the moving average filter is about 0.5 at
% this frequency.

%unit sample response and frequency response
h = ones(5,5)/25;
[H,fx,fy] = freqz2(h,64,64,[1 1]);
[Fx,Fy] = meshgrid(fx,fy);
%contour plot
figure;
v=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.2
[C,h1]=contour(Fx,Fy,abs(H),v); % generate the contour plot, including values to label contours
axis square
v1 = 0:.2:1.;
clabel(C,h1,v1,'FontSize',6);           %label the contours  
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 
%set(gca,'XAxisLocation','top');      %Xaxis labels on top
colormap([0 0 0]);     % use black only
hold on
% plot a circle of radius 1/8
rectangle('Position',[-0.125, -0.125, 0.25, 0.25],'Curvature',[1,1],'EdgeColor','r','LineWidth',1);
axis square
set(gca,'ydir','reverse'); 
%set(gca,'XAxisLocation','top');      %Xaxis labels on top
colormap([0 0 0]);     % use black only
hold off
set(gcf,'Color',[1 1 1]);
%
%perspective plot
figure;
mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(H(1:2:64,1:2:64)));            % generate the perspective plot
colormap([0 0 0]);     % use black only
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 
%write the plot to a file
set(gcf,'Color',[1 1 1]);
%
%Filter the barbara image
A = im2double(imread('barb512_gray.tif'));
B = conv2(A,h,'same');
figure;
imshow(B)

% Gaussian filter
% To meet the specs, we need r = 1.06 and c = 0.1416 (see notes). 
L=6;
x = -L:1:L;
[X,Y]= meshgrid(x);
r = 1.06;
c = .1416;
hg = c*exp(-(X.^2 + Y.^2)/(2*r^2));
[HG,fx,fy] = freqz2(hg,64,64,[1 1]);
[Fx,Fy] = meshgrid(fx,fy);
v=0:.1:1.;             % contours will be from 0 to 1 in steps of 0.1
figure;
[C,h1]=contour(Fx,Fy,abs(HG),v); % generate the contour plot, including values to label contours
axis square
v1 = 0:.2:1.;
clabel(C,h1,v1,'FontSize',6);           %label the contours  
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 
%set(gca,'XAxisLocation','top');      %Xaxis labels on top
colormap([0 0 0]);     % use black only
hold on
% plot a circle of radius 1/8
rectangle('Position',[-0.125, -0.125, 0.25, 0.25],'Curvature',[1,1],'EdgeColor','r','LineWidth',1);
axis square
set(gca,'ydir','reverse'); 
%set(gca,'XAxisLocation','top');      %Xaxis labels on top
colormap([0 0 0]);     % use black only
hold off
set(gcf,'Color',[1 1 1]);
%
%perspective plot
figure;
mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(HG(1:2:64,1:2:64)))            % generate the perspective plot
colormap([0 0 0]);     % use black only
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse');

set(gcf,'Color',[1 1 1]);
%
%Filter the barbara image
D = conv2(A,hg,'same');
figure,imshow(D)