clear all;
close all;
%---------------------------------------------------------------
%Initialize colormap
mcomp = 0 : 1.0 ./ 255 : 1.0; % output colors lie in [0,1]
map = [mcomp' mcomp' mcomp'];

M = 256;  % size of image, must be even
C = 8; % sampling interval in x and y, must be even
gamma = 1.0;
gamma_inv = 1/gamma;
index = (0:1:M-1);
y = ones(M,1)*index;
x = y';
%degrees = 10;
%theta = pi*degrees/180;
%P = 10;
%rho = 1/P;
%u = rho*cos(theta)
%v = rho*sin(theta)
u = 1/32
v = 1/8
g = 0.5 + 0.5*cos(2*pi*(u*x + v*y));
figure(1);
subplot(2,3,1)
image(g.^gamma_inv*255);
colormap(map);
axis('off');
axis('image');
title('(a)')

g = fftshift(g);
G = fft2(g);
G = fftshift(G);
Gmag = abs(G);
Gmag = conv2(Gmag, ones(2,2), 'same');
Gmax = max(max(G));
figure(1);
subplot(2,3,4)
image(Gmag*255/Gmax);
colormap(map);
axis('off');
axis('image');
title('(d)')
clear G;
clear Gmag;
%----------------------------------------
%Downsample, and then upsample by (C,D)
g_s = zeros(M,M);
g_s(1:C:M,1:C:M) = g(1:C:M,1:C:M);
g_s_disp = conv2(g_s, ones(2,2), 'same');
figure(1);
subplot(2,3,2)
image(g_s_disp*255);
colormap(map);
axis('off');
axis('image');
title('(b)')

g_s = fftshift(g_s);
G_s = fft2(g_s);
G_s = fftshift(G_s);
G_smag = abs(G_s);
G_smag = conv2(G_smag, ones(2,2), 'same');
G_smax = max(max(G_s));

%-----------------------------------------
% Show baseband region
U = M/C;
u0 = M/2 - U/2;
v0 = M/2 - U/2;
u1 = u0+U-1;
v1 = v0;
u2 = u1;
v2 = v0+U-1;
u3 = u0;
v3 = v2;
for k = 0:U-2,	
	G_smag(u0+k+1,v0+1) = G_smax;
	G_smag(u1+1,v1+k+1) = G_smax;
	G_smag(u2-k+1,v2+1) = G_smax;
	G_smag(u3+1,v3-k+1) = G_smax;
end

figure(1);
subplot(2,3,5)
image(G_smag*255/G_smax);
colormap(map);
axis('off');
axis('image');
title('(e)')

clear g_s;
clear G_smag;
%-----------------------------------------
% Reconstruct
G_r = zeros(M,M);
G_r(M/2-U/2+1:M/2+U/2,M/2-U/2+1:M/2+U/2) = ...
C.*C.*G_s(M/2-U/2+1:M/2+U/2,M/2-U/2+1:M/2+U/2);
G_rmag = abs(G_r);
G_rmag = conv2(G_rmag, ones(2,2), 'same');
G_rmax = max(max(G_r));

figure(1);
subplot(2,3,6)
image(G_rmag*255/G_rmax);
colormap(map);
axis('off');
axis('image');
title('(f)')

clear G_s;
clear G_rmag;

G_r = fftshift(G_r);
g_r = ifft2(G_r);
g_r = fftshift(g_r);

figure(1);
subplot(2,3,3)
image(g_r.^gamma_inv*255);
colormap(map);
axis('off');
axis('image');
title('(c)')