aliasing2.m

aliasing in matlab by lrmaster,aliasing in matlab by lrmaster

clear all;
close all;

%%%%%%%%%%%%%
% Make up a digital camera
%
% define a focal point
% define a viewing plane
% define a window of pixels
focus = [0 0 -1]'; % define focal point
vwcnt = [0 0 0]'; % centre of view plane
vwup = [0 1 0]';  % view plane up direction
vwrt = [1 0 0]';   % view plane right direction
vwfd = [0 0 1]';   % view plane forward direction - normal to plane
Npix = 384;          % the number of pixels
window = zeros( Npix, Npix, 3 );  % the rgb pixels
[pixx,pixy] = meshgrid( 0:(Npix-1) ); % define the location of pixel centres on the viewplane
pixx = pixx - Npix/2 + 0.5; % shift the centre to a sensible place...
pixy = pixy - Npix/2 + 0.5; % ...so the optical axis runs right trhough the midel of the view plane
pixy = flipud(pixy); % make the y values the right way around
% scale all values do the window is of unit size
pixx = pixx / (2*max(max(abs(pixx))));
pixy = pixy / (2*max(max(abs(pixy))));

%%%%%%%%%%%%%%
% Make up a chessboard floor
%
floorcentre = [0 -10 0]';
floornormal = [0 1 0]';
floorup = [0 0 1]';
floorrt = [1 0 0]';

% scale and rotate the pattern
a = pi/3;
R = [ cos(a) 0 sin(a);
  0 1 0;
  -sin(a) 0 cos(a)
];
S = [0.5 0 0; 0 1 0; 0 0 0.5];
M = R*S;
floorup = M * floorup;
floorrt = M * floorrt;

for j = 1:Npix
  for i = 1:Npix
    % get position of the pixel in the eye frame
    x = vwcnt + pixx(i,j) * vwrt + pixy(i,j) * vwup;
    % get the ray direction
    ray = x - focus;
    % get the cosine of the angle or floor normal to ray
    s = ray'*floornormal;
    if s < eps
      % find the distance intersection of the ray with the floor
      s = ((floorcentre - focus)' * floornormal) / s;
      if s > 0
        % now compute the actual intersection
        y = focus + s*ray;
        % the y is in the eye frame, move it into the floor frame
        y = y - floorcentre; % shift to floor origin
        y = [floorrt floorup]' * y; % compute coordinates in floor frame
        % now decide the area that the point is in
        u = rem(y,2);
        u = rem( u+2, 2 );
        u = u > 1;
        switch( u(1) +  2*u(2) )
        case 0
          window(i,j,1) = 0;
          window(i,j,2) = 1;
          window(i,j,3) = 0;
        case 1
          window(i,j,1) = 0;
          window(i,j,2) = 0.5;
          window(i,j,3) = 0;
        case 2
          window(i,j,1) = 0;
          window(i,j,2) = 0.5;
          window(i,j,3) = 0;
        case 3
          window(i,j,1) = 1;
          window(i,j,2) = 1;
          window(i,j,3) = 1;
        end
      else
        % the intersection is behind the focal point, it cannot be seen
      end
    else
      % the ray is parallel to the floor
    end
  end
end

%window(:,:,1) = pictNorm( window(:,:,1) );
%window(:,:,2) = pictNorm( window(:,:,2) );

imshow( window );