% fig_sample_recon.m
% illustrate sampling and then reconstruction (interpolation)

if ~isvar('x1') || 0 % synthetic
	ig = image_geom('nx', 2^10, 'fov', 24);
	x1 = ellipse_im(ig, 'shepplogan-emis', 'oversample', 3);
	clf, im(ig.x, ig.y, x1)
	nx = ig.nx;
	ny = ig.ny;
end

if ~isvar('x1') || 0
	file = '2016_06_23_hallway.jpg';
	x0 = imread(file);
	image(x0)
	x1 = rgb2gray(x0);
%	x1 = x1([1:1536]+150,[1:2560]+100);
	x1 = x1([1:1536]+150,[1:2048]+400);
	x1 = flipdim(x1,1);
	im(x1)
	pr factor(size(x1,1))
	pr factor(size(x1,2))
	x1 = downsample2(x1, 2);
	[nx ny] = size(x1);
end

du = 1/(ig.nx*ig.dx);
dv = 1/(ig.ny*ig.dy);
ux = [-nx/2:nx/2-1]/nx * du;
uy = [-ny/2:ny/2-1]/ny * dv;

%ox = [-nx/2:nx/2-1]/nx * 2 * pi;
%oy = [-ny/2:ny/2-1]/ny * 2 * pi;
fun = @(x) max(20*log10(abs(x) / max(abs(x(:)))), -140);

down = 4;
x2 = x1(1:down:nx,1:down:ny);
x2i = zeros(nx,ny);
x2i(1:down:nx,1:down:ny) = x2;

ir_fontsize title 10
ir_fontsize im_axes 8
ir_fontsize label 11

%clim = [0 255];
clim = [0 10]; % shepp

ulim = [-1 1] * 0.02;
utic = [-1 0 1] * 0.02;

%imd = @(n, xx) im(5, 'notick', x, y, xx(ix,iy), clim);
%imd = @(n, xx) im(5, ig.x(ix), ig.y(iy), xx(ix,iy), clim);
imd = @(n, xx) im(n, ig.x, ig.y, xx, clim);
imdn = @(n, xx) im(n, 'notick', ig.x, ig.y, xx, clim);

%imf = @(n, ff) im(n, ox, oy, fun(ff));
imf = @(n, ff) im(n, ux, uy, fun(ff));
imfn = @(n, ff) im(n, 'notick', ux, uy, fun(ff));

if 1
	box = zeros(nx,ny);
	box(nx/2+1+[-nx/(2*down):nx/(2*down)],ny/2+1+[-ny/(2*down):ny/(2*down)]) = 1;
	X2r = box .* fftshift(fft2(x2i));

	x2r = down^2 * ifft2(ifftshift(X2r));

	X1f = box .* fftshift(fft2(x1));
	x1f = ifft2(ifftshift(X1f));

	im plc 3 4

	im(1, ig.x, ig.y, x1, clim)
	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
	labels('$x$ [cm]', '$y$')
	titlef 'Original image \gxy'

	im(2, 'notick', ig.x, ig.y, x2i, clim)
%	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
%	labels('$x$', '$y$')
	titlef 'Image samples $g_s(x,y)$'

	im(3, 'notick', ig.x, ig.y, x2r, clim)
%	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
%	labels('$x$', '$y$')
	titlef 'Recovered image'

	im(4, 'notick', ig.x, ig.y, x1f, clim)
%	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
%	labels('$x$', '$y$')
	titlef 'Prefiltered image'

	ix = 350+[-50:50];
	iy = 500+[-50:50];
	ix = 480+[-50:50];
	iy = 830+[-50:50];

	imd(5, x1)
	axisx([-1.9 .4]), axisy([-8.4 -6.2])
	xtick([-1.9 .4]), ytick([-8.4 -6.2])
	labels('$x$ [cm]', '$y$')
	titlef '(zoom)'

	imdn(6, x2i)
	axisx([-1.9 .4]), axisy([-8.4 -6.2])
%	xtick([-1.9 .4]), ytick([-8.4 -6.2])
%	labels('$x$', '$y$')
	titlef '(zoom)'

	imdn(7, x2r)
	xlim([-1.9 .4]), ylim([-8.4 -6.2])
%	xtick([-1.9 .4]), ytick([-8.4 -6.2])
%	labels('$x$', '$y$')
	titlef '(zoom)'

	imdn(8, x1f)
	xlim([-1.9 .4]), ylim([-8.4 -6.2])
%	xtick([-1.9 .4]), ytick([-8.4 -6.2])
%	labels('$x$', '$y$')
	titlef '(zoom)'


	imf(9, fftshift(fft2(x1)))
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim), xtick(utic), ytick(utic)
	xlabelf '\kx [1/cm]', ylabelf '\ky',
	titlef 'Original image spectrum'

	imfn(10, fftshift(fft2(x2i)))
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim)
%	xtick(utic), ytick(utic)
%	xlabelf '\kx', ylabelf '\ky',
	titlef 'Sampled image spectrum'

	imfn(11, X2r)
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim)
%	xtick(utic), ytick(utic)
%	xlabelf '\kx', ylabelf '\ky',
	titlef 'Recovered image spectrum'

	imfn(12, X1f)
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim)
%	xtick(utic), ytick(utic)
%	xlabelf '\kx', ylabelf '\ky',
	titlef 'Prefiltered image spectrum'

%	ir_savefig fig_sample_recon1
return
end

if 1
	x2s = upsample_rep(x2, down);
	im plc 3 4

	x2p = zeros([nx ny]);
	x2p(480,830) = 9;
	x2p = conv2(x2p, ones(down), 'same');

	im(1, ig.x, ig.y, x1, clim)
	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
	labels('$x$ [cm]', '$y$')
	titlef 'Original image \gxy'

	im(2, 'notick', ig.x, ig.y, x2i, clim)
%	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
%	labels('$x$', '$y$')
	titlef 'Image samples $g_s(x,y)$'

	im(4, 'notick', ig.x, ig.y, x2s, clim)
%	xtick([-1 1]*ig.fov/2), ytick([-1 1]*ig.fov/2)
%	labels('$x$', '$y$')
	titlef 'Interpolated image'


	imd(5, x1)
	axisx([-1.9 .4]), axisy([-8.4 -6.2])
	xtick([-1.9 .4]), ytick([-8.4 -6.2])
	labels('$x$ [cm]', '$y$')
	titlef '(zoom)'

	imdn(6, x2i)
	axisx([-1.9 .4]), axisy([-8.4 -6.2])
%	xtick([-1.9 .4]), ytick([-8.4 -6.2])
%	labels('$x$', '$y$')
	titlef '(zoom)'

	imdn(7, x2p)
	axisx([-1.9 .4]), axisy([-8.4 -6.2])
%	xtick([-1.9 .4]), ytick([-8.4 -6.2])
%	labels('$x$', '$y$')
	titlef 'Interp PSF'

	imdn(8, x2s)
	axisx([-1.9 .4]), axisy([-8.4 -6.2])
%	xtick([-1.9 .4]), ytick([-8.4 -6.2])
%	labels('$x$', '$y$')
	titlef '(zoom)'


	imf(9, fftshift(fft2(x1)))
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim), xtick(utic), ytick(utic)
	xlabelf '\kx [1/cm]', ylabelf '\ky',
	titlef 'Original image spectrum'

	imfn(10, fftshift(fft2(x2i)))
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim)
%	xtick(utic), ytick(utic)
%	xlabelf '\kx', ylabelf '\ky',
	titlef 'Sampled image spectrum'

	imfn(11, fftshift(fft2(x2p)))
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim)
%	xtick(utic), ytick(utic)
%	xlabelf '\kx', ylabelf '\ky',
	titlef 'Interpolator Freq. Resp.'

	imfn(12, fftshift(fft2(x2s)))
	%axis square, axis_pipi
	xlim(ulim), ylim(ulim)
%	xtick(utic), ytick(utic)
%	xlabelf '\kx', ylabelf '\ky',
	titlef 'Interpolated image spectrum'

%	ir_savefig fig_sample_recon2
end

function labels(x, y)
xlabelf(x)
axisxnudge(10)
ylabelf(y)
axisynudge(-12)
end
