% projections.m, 3D to 2D
% run this program several times, random data.
clear

% Generate the space on which to project
v1=randn(3,1); v2=randn(3,1); A=[v1 v2];

% Construct the projection matrix and check for key properties
P=A*inv((A'*A))*A';
norm(P-P','fro'), norm(P*P-P,'fro')

% Create an ON-basis for the plane and generate plane-plot
[q,r]=qr(A,0);
step=0.2; width=3; index=1;
for l=-width:step:width
	p1=l*q(:,1)-width*q(:,2);
	p2=l*q(:,1)+width*q(:,2);
	p3=(l+step/2)*q(:,1)+width*q(:,2);
	p4=(l+step/2)*q(:,1)-width*q(:,2);
	planeplot(:,index)=p1;
	planeplot(:,index+1)=p2;
	planeplot(:,index+2)=p3;
	planeplot(:,index+3)=p4;
	index=index+4;
end
ax=[ -width width -width width -width width];

% Generate some random vectors and project them
M=50;
dum=randn(3,M);
dummer=P*dum;

% plot some examples
figure(1), clf, hold on, axis(ax), axis equal, view(3)
plot3([0 v1(1)], [0 v1(2)], [0 v1(3)], 'b')
plot3([0 v2(1)], [0 v2(2)], [0 v2(3)], 'b')
for l=1:5
	plot3(dum(1,l),dum(2,l),dum(3,l),'rp')
	plot3(dummer(1,l),dummer(2,l),dummer(3,l),'gp')
	plot3([dum(1,l) dummer(1,l)], [dum(2,l) dummer(2,l)],...
		 [dum(3,l) dummer(3,l)], 'y')
end
plot3(planeplot(1,:), planeplot(2,:), planeplot(3,:), 'k')
title('Five points and their projections')

% plot all projected points
figure(2), clf, hold on, axis(ax), axis equal, view(3)
plot3([0 v1(1)], [0 v1(2)], [0 v1(3)], 'b')
plot3([0 v2(1)], [0 v2(2)], [0 v2(3)], 'b')
for l=1:M
	plot3(dummer(1,l),dummer(2,l),dummer(3,l),'gp')
end
plot3(planeplot(1,:), planeplot(2,:), planeplot(3,:), 'k')
title('Fifty projected points')