% Chapter 11, Fig.11.10
% Eric Dubois, updated 2019-01-14
% Downsample Barbara image by two in each dimension
% Images and spectral estimates 
% Required data: barb512_gray.tif, antialias_qs_coeff.mat
% Required functions: WelchPDS_est
% antialias_qs_coeff.mat is generated using mdsp_web_changesamp_design_qs_antialiasing_filter.m
clear all, close all;
N=64; % number of points of spectral estimate
A = im2double(imread('barb512_gray.tif'));

%power density spectrum of the original Barbara image, in dB
BarbPDS = WelchPDS_est(A,N);
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
BarbPDS_in_dB = 20*log10(BarbPDS);
ground = -100;
BarbPDS_in_dB(BarbPDS_in_dB<ground) = ground;
maxplot = ceil(max(max(BarbPDS_in_dB))/10)*10;
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),BarbPDS_in_dB(1:2:N,1:2:N));
xlabel('u (c/ipx)'), ylabel('v (c/ipx)')
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]);

%set three quarters of samples to zero to get spectrum of subsampled image.
%and plot power density spectrum

B_zero_ins = A;
B_zero_ins(2:2:size(A,1),2:2:size(A,2)) =0;
figure, imshow(B_zero_ins);

BarbDSPDS = WelchPDS_est(B_zero_ins,N);
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
BarbDSPDS_in_dB = 20*log10(BarbDSPDS);
ground = -100;
BarbDSPDS_in_dB(BarbDSPDS_in_dB<ground) = ground;
maxplot = ceil(max(max(BarbDSPDS_in_dB))/10)*10;
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),BarbDSPDS_in_dB(1:2:N,1:2:N));
xlabel('u (c/ipx)'), ylabel('v (c/ipx)')
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]);


%Downsampled image with no prefilter and save
B = A(1:2:size(A,1),1:2:size(A,2));
figure, imshow(B);

%Prefilter the Barbara image with a prefilter designed with
%mdsp_fig_changesamp_design_qs_antialiasing_filter.m
load('antialias_qs_coeff.mat');
A_aapf = imfilter(A,h_aa,'same','symmetric');
figure, imshow(A_aapf);
%apply zero insertion and compute and plot spectra before and after
%prefiltering
B_aapf_zero_ins = A_aapf;
B_aapf_zero_ins(2:2:size(A,1),2:2:size(A,2)) =0;
figure, imshow(B_aapf_zero_ins);

%power density spectrum of the prefiltered Barbara image, in dB
BarbPFPDS = WelchPDS_est(A_aapf,N);
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
BarbPFPDS_in_dB = 20*log10(BarbPFPDS);
ground = -100;
BarbPFPDS_in_dB(BarbPFPDS_in_dB<ground) = ground;
maxplot = ceil(max(max(BarbPFPDS_in_dB))/10)*10;
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),BarbPFPDS_in_dB(1:2:N,1:2:N));
xlabel('u (c/ipx)'), ylabel('v (c/ipx)')
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]);

BarbPFDSPDS = WelchPDS_est(B_aapf_zero_ins,N);
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
BarbPFDSPDS_in_dB = 20*log10(BarbPFDSPDS);
ground = -100;
BarbPFDSPDS_in_dB(BarbPFDSPDS_in_dB<ground) = ground;
maxplot = ceil(max(max(BarbPFDSPDS_in_dB))/10)*10;
figure;
mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),BarbPFDSPDS_in_dB(1:2:N,1:2:N));
xlabel('u (c/ipx)'), ylabel('v (c/ipx)')
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]);

%Downsample the prefiltered image
B_aapf = A_aapf(1:2:size(A,1),1:2:size(A,2));
figure, imshow(B_aapf);

C = imresize(A,.5,'bicubic');
figure, imshow(C);