flow, random points on a plane

% calculate the scene - a field of points
clear
close all
N = 2000;
x = rand(N,1)*40-20;
y = rand(N,1)*80-20;
z = zeros(size(x));
% draw the scene
figure('Position',[320 240 320 240]);
plot3(x,y,z,'ko');
axis equal
axis off
camproj('perspective');
campos([0 -20 5]);
camtarget([0 -10 0]);
camva(30);
pos = campos();
tgt = camtarget();
% move the camera
delta = [0 0.2 0];
aviobj = avifile('rand_flow1.avi','fps',20,'compression','Cinepak');
for k=1:81
    campos(pos+(k-1)*delta);
    camtarget(tgt+(k-1)*delta);
    drawnow;
    F = getframe(gcf);
    aviobj = addframe(aviobj,F);
end
aviobj = close(aviobj);
% do flow without MATLAB 3D graphpics


figure(2)
up = [0 0 1];
vn = pos - tgt;
vz = vn/norm(vn);
vx = cross(up,vz);
vx = vx/norm(vx);
vy = cross(vz,vx);
R = [vx; vy; -vz;];
T2 = [R  -R*pos'];

% test
disp('transform target point:');
disp(T2*[tgt 1]')

xp = x(:);
yp = y(:);
zp = z(:);

xform = T2*[xp'; yp'; zp'; ones(size(zp'))];
va = 30;
scl = 1/tan(0.5*va*pi/180);
idx = find(xform(3,:)>0); % select points in front of camera
xpts = xform(1,idx)./xform(3,idx)*scl;
ypts = xform(2,idx)./xform(3,idx)*scl;

plot(xpts,ypts,'ko');
axis([-1 1 -1 1]);
transform target point:
         0
         0
   11.1803

% traditional optical flow diagram
yp = y(:) - 0.2;
xform = T2*[xp'; yp'; zp'; ones(size(zp'))];
xnew = xform(1,idx)./xform(3,idx)*scl;
ynew = xform(2,idx)./xform(3,idx)*scl;
quiver(xpts,ypts,xnew-xpts,ynew-ypts,10,'k');
axis([-1 1 -1 1]);