% ------------------------------------------------------------------------
%   ECE 178 - Tutorial 3
%   by Marco Zuliani - 01/31/2002
%
%   DESCRIPTION
%   In this example we will explore the effects of the aliasing phenomenon
%   in a two dimensional case
% ------------------------------------------------------------------------


% let's clear the memory of the workspace
clear all
% let's clear the command window
clc
% and let's close all the figures
close all


% let's define the spatial base for the pseudo-continuous signal (note
% how the time increment between samples is very small!
x_min = -0.5;
y_min = -0.5;
x_max = 0.5;
y_max = 0.5;
x = x_min:2*1e-3:x_max;
y = x_min:2*1e-3:y_max;
[X Y] = meshgrid(x, y);
u = 4;
v = 5;
I = sin(2*pi*(u*X+v*Y));
clear X Y


% and visualize the result
figure
subplot(1, 2, 1)
imagesc(I);
colormap(gray);
axis equal
xlabel('x')
ylabel('y');
title('Original image');
drawnow


% let's sample the image
us = 13;
vs = 15;
deltax = 1/us;
deltay = 1/vs;
xs = x_min:deltax:x_max;
ys = y_min:deltay:y_max;
[Xs Ys] = meshgrid(xs, ys);
Is = sin(2*pi*(u*Xs+v*Ys));
clear Xs Ys


% and now let's try to reconstruct it
Ir = zeros(size(I));
Nx = size(Is, 2);
Ny = size(Is, 1);
for n = 1:Nx
    for m = 1:Ny;
        
        % generate sinc x
        den = pi*(x - (n-0.5*Nx)*deltax)/deltax;
        num = sin(den);
        % take care of zero denominators (sinc(0) is redifined to be 1)
        ind = find(den == 0);
        if ~isempty(ind)
            den(ind) = 1;
            num(ind) = 1;
        end;
        sincx = num ./ den; 
        sincx = repmat(sincx, [size(I, 1), 1]);
        
        % generate sinc y
        den = pi*(y - (m-0.5*Ny)*deltay)/deltay;
        num = sin(den);
        % take care of zero denominators (sinc(0) is redifined to be 1)
        ind = find(den == 0);
        if ~isempty(ind)
            den(ind) = 1;
            num(ind) = 1;
        end;
        sincy = num ./ den;
        sincy = repmat(sincy', [1, size(I, 2)]);
        
        Ir = Ir + Is(m, n) .* sincx .* sincy;
        
    end
end


% visualize the reconstructed image
subplot(1, 2, 2)
imagesc(Ir);
colormap(gray);
axis equal
xlabel('x')
ylabel('y');
title('Reconstructed image');