I was looking for an easy way to derive a projection matrix corresponding to the camera of images rendered with blender. Given such a 3×4 projection matrix P, I want to define a 3D point in the blender scene, multiply it with the matrix and paint it into the existing image. For example, using Matlab to draw the scene origin as a red plus sign onto the image, I want to write
imshow('blender_img.jpg');
hold on;
X = [0 0 0 1]';
x = P*X;
x = x(1:2)/x(3);
plot(x(1),x(2),'r+');
hold off;
Deriving P turns out to be not documented at all, and besides knowledge on the computer graphics pipelines requires some trial and error. Today I managed to solve the problem, so here's the solution.

First of all, collect the camera position and orientation from the "Transform Properties" dialog, as shown in the image above. In Matlab you get the translation vector and rotation matrix via
X0 = [2 -10 4]';
o = 70*pi/180;
p = 15*pi/180;
k = -10*pi/180;
Ro=[1 0 0; 0 cos(o) -sin(o); 0 sin(o) cos(o)];
Rp=[cos(p) 0 sin(p); 0 1 0; -sin(p) 0 cos(p)];
Rk=[cos(k) -sin(k) 0; sin(k) cos(k) 0; 0 0 1];
R = Rk*Rp*Ro;
Now we need to derive the calibration matrix. Unfortunately, the focal length value in blender is neither documented nor standard. However, it is possible to display the lens angle alpha in degrees, referring to the wider image dimension, as shown in the image below. In the example we have alpha=49.13.

From the image resolution, which is here assumed to be 800×600, we get the focal length in pixel from fl=-400/(49.13*pi/180/2). Observe the minus!
The final projection matrix P is obtained in Matlab using
fl = -400/(49.13*pi/180/2);
K = [fl 0 400; 0 fl 300; 0 0 1];
P = K*diag([1 -1 1])*R'*[eye(3) -X0];
Note the flipping of the y axis which results in a left handed system. A
student is planning to write an export script for blender to write
P to file, which I'll post here.