%Chapter 10, Example 10.5, Fig 10.18, 10.19
%Eric Dubois, updated 2018-12-21
% Low pass filter with frequency domain constraints to filter halftone image
% Use a zero phase filter. Do not use quadrantal symmetry because the spec
% is skewed.
% Result without constraints can also be generated by using filter_design_zp_unc
% Required functions: ideal_response_zp, filter_design_zp_unc,
%    filter_design_zp_con, WelchPDS_est
% Required data files: halftone_lpf_specs_40_60.txt,
%    halftone_lpf_constraints.txt, eric_prize_1963_crop_gray.jpg
close all, clear all;
%Generate ideal filter response
Npt=64;
[f1,f2] = freqspace(Npt);
[F1,F2] = meshgrid(.5*f1,.5*f2);
Hd = ideal_response_zp('halftone_lpf_specs_40_60.txt',Npt);
% perspective view of the ideal response
%subsample the response to get a suitable plot
figure;
mesh(F1(1:2:Npt,1:2:Npt),F2(1:2:Npt,1:2:Npt),abs(Hd(1:2:Npt,1:2:Npt)));            % generate the perspective plot
view(-18,40);
axis([-.5, .5, -.5, .5, 0, 1.2]);
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]);
%Design filter
%h = filter_design_zp_con(Hd,21,21,'halftone_lpf_constraints.txt');
h = filter_design_zp_unc(Hd,21,21);
%Compute and plot the frequency response
[H,fx,fy] = freqz2(h,64,64,[1 1]);
[Fx,Fy] = meshgrid(fx,fy);
%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
view(-18,40);
colormap([0 0 0]);     % use black only
axis([-.5 .5 -.5 .5 0 1.2]);
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 

set(gcf,'Color',[1 1 1]);

%Filter the halftone image with the lowpass filter
INPUT = im2double(imread('eric_prize_1963_crop_gray.jpg'));
OUTPUT = imfilter(INPUT,h,'same','symmetric');
figure, imshow(INPUT);
figure, imshow(OUTPUT);

N=64;
spect_out = WelchPDS_est(OUTPUT,N);
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
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
spect_out_dB = 20*log10(spect_out);
ground = -80;
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]);