function [alpha, beta, theta, u0, v0, R, t] = calibrate(f2D,f3D)
%% 
%% computes the intrinsic (camera) and extrinsic (world) parameters
%% given real-world coordinates in homogenous coords (f3D) of some
%% item in space, and corresponding pixel locations in f2d
%%

M = projM(f2D,f3D) % first estimate the projection matrix

% M = (A b)
A = M(:,[1:3]);
b = M(:,4);
% then recover the parameters
a1 = A(1,:)'; a2 = A(2,:)'; a3 = A(3,:)';

% eqtns 3.11 in forsythe + ponce
rho = 1/sqrt(sum(a3.^2))
r3 = rho * a3
u0 = rho.^2 * (a1' * a3)
v0 = rho.^2 * (a2' * a3);

% equations 3.13 in forsythe + ponce
theta = acos( -(cross(a1, a3)'*cross(a2,a3)) / ( sqrt(sum(cross(a1,a3).^2)) * sqrt(sum(cross(a2,a3).^2)) ) );
alpha = rho.^2 * sqrt(sum(cross(a1,a3).^2))*sin(theta);
beta = rho.^2 *  sqrt(sum(cross(a2,a3).^2))*sin(theta);

% eqtns 3.14 in forsythe + ponce
%r1 = 1 /(norm(cross(a2,a3))) * cross(a2,a3);
r1 = rho ^2 * sin(theta) / beta * cross(a2,a3);
r2 = cross(r3,r1);

R = [r1'; r2'; r3']
K = [alpha -alpha*cot(theta) u0; 0 beta/sin(theta) v0; 0 0 1]; %rho * A * inv(R);
t = rho * inv(K)*b;