%Chapter 10, Example 10.3, Fig. 10.8, 10.9, 10.11, 10.12
%Eric Dubois, updated 2018-12-20
% Multi bandstop filter to attenuate halftone pattern
% Assume the lattice is Z^2. 
% Do a spectral analysis of the halftone picture to determine the main
% frequencies of the halftone pattern
%Required functions: WelchPDS_est

% Read the input image
close all; clear all;
INPUT = im2double(imread('eric_prize_1963_crop_gray.jpg'));
%Baseband Gaussian filter
L = 8;  %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
r = 3.0; % parameter of the Gaussian filter, determined empirically
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
hLP = h1*c;
%Enter the frequencies of the halftone interference within a unit cell of 
%the lattice Z^2 in c/px
FB(:,1) = [0.14737 0.16140];
FB(:,2) = [0.16140, - 0.14737];
FB(:,3) = FB(:,1) + FB(:,2);
FB(:,4) = FB(:,1) - FB(:,2);
FB(:,5) = 2*FB(:,1);
FB(:,6) = 2*FB(:,2);
FB(:,7) = FB(:,1) + FB(:,3);
FB(:,8) = FB(:,1) + FB(:,4);
FB(:,9) = FB(:,2) + FB(:,3);
FB(:,10) = FB(:,2) - FB(:,4);
%Generate the multiband bandpass filter; use the first NF frequencies,
% where NF can be 4, 6, or 10
NF = 10;
hBP = 0*hLP;
for K = 1:4;
    hBP = hBP +  2 * hLP .* cos(2*pi*(X*FB(1,K)+Y*FB(2,K)));
end
if NF > 4
    for K = 5:6;
      hBP = hBP +  2 * hLP .* cos(2*pi*(X*FB(1,K)+Y*FB(2,K)));
    end
end
if NF > 6
    % find the extra gain due to the overlap of nearby components near the
    % edge of the unit cell
    hBP_7_9 = 2 * hLP .* (cos(2*pi*(X*FB(1,7)+Y*FB(2,7))) + ...
        cos(2*pi*(X*FB(1,9)+Y*FB(2,9))));
    HBP_7_9 = freqz2(hBP_7_9);
    maxgain_7_9 = max(max(HBP_7_9));
    for K = 7:10
        hBP = hBP +  (2/maxgain_7_9) * hLP .* cos(2*pi*(X*FB(1,K)+Y*FB(2,K)));
    end
end
%generate the bandstop filter
hBS = zeros(2*L+1,2*L+1);
hBS(L+1,L+1)=1;
hBS = hBS - hBP;
%Filter the halftone image with the bandstop filter
OUTPUT = imfilter(INPUT,hBS,'same','symmetric');
figure, imshow(INPUT);
figure, imshow(OUTPUT);
%Plot spectra of input and output images
% Size the mesh plots to fit side by side
N=64;
spect_in = WelchPDS_est(INPUT,N);
spect_out = WelchPDS_est(OUTPUT,N);
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
%Linear scale
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),spect_in(1:2:N,1:2:N));
xlabel('u (c/px)'), ylabel('v (c/px)')
zlabel('Power Density')
colormap([0 0 0]);     % use black only
set(gca,'ydir','reverse');  %make v axis point downward
set(gcf,'Color',[1 1 1]);
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),spect_out(1:2:N,1:2:N));
xlabel('u (c/px)'), ylabel('v (c/px)')
zlabel('Power Density')
colormap([0 0 0]);     % use black only
set(gca,'ydir','reverse');  %make v axis point downward
set(gcf,'Color',[1 1 1]);
%dB scale
spect_in_dB = 20*log10(spect_in);
ground = -60;
spect_in_dB(find(spect_in_dB<ground)) = ground;
maxplot = ceil(max(max(spect_in_dB))/10)*10;
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),spect_in_dB(1:2:N,1:2:N));
xlabel('u (c/px)'), ylabel('v (c/px)')
zlabel('Power Density (dB)')
colormap([0 0 0]);     % use black only
set(gca,'ydir','reverse');  %make v axis point downward
set(gcf,'Color',[1 1 1]);
spect_out_dB = 20*log10(spect_out);
ground = -60;
spect_out_dB(find(spect_out_dB<ground)) = ground;
maxplot = ceil(max(max(spect_out_dB))/10)*10;
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),spect_out_dB(1:2:N,1:2:N));
xlabel('u (c/px)'), ylabel('v (c/px)')
zlabel('Power Density (dB)')
colormap([0 0 0]);     % use black only
set(gca,'ydir','reverse');  %make v axis point downward
set(gcf,'Color',[1 1 1]);
% Contour plot of input spectrum in dB
figure;
vlev=ground:10:maxplot;             % contours will be in steps of 10 dB
[Con,hval]=contour(u,v,spect_in_dB,vlev); % generate the contour plot, including values to label contours
axis square;  %make the plot square
xlabel('Horizontal frequency (c/px)'), ylabel('Vertical frequency (c/px)')
clabel(Con,hval,'FontSize',6)             %label the contours
set(gca,'ydir','reverse');  %make v axis point downward
colormap([0 0 0]);     % use black only
%write the plot to a file
set(gcf,'Color',[1 1 1]);