Calculating Projection Matrices for MATLAB/Simulink Virtual Reality Toolbox

I was a bit stumped this week, trying to calculate the projection matrices which MATLAB/Simulink uses when rendering 3D worlds (X3D/VRML). I will document it here, perhaps it will be helpful to someone. I don't know anything about computer graphics, so I don't know the correct terminology for this sort of thing. Let's start with a simple bit of MATLAB.

campos = [ 0 0 300 ]; % camera position
rendW = 384;          % output width
rendH = 384;          % output height

world = vrworld('path/to/vrml/file.wrl');
open(world)
canvas = vrfigure(world);
set(canvas, ...
       'CameraPosition', campos, ...
       'Position', [ 0 0 rendW rendH ])
I = capture(canvas);
close(canvas)
close(world)

We can extract the vertices of the object in a similar way to below, but be careful of the field name, since it will change depending on the world file. While the world is open:

P = getfield(vrnode(world, 'di3d-COORD'), 'point');

So, these points have had no transformations applied to them, they are not projected, they are as they would be in the VRML/X3D file. So, let's begin by calculating our camera matrix.

C = [ 1 0 0 -campos(1)
      0 1 0 -campos(2)
      0 0 1 -campos(3)
      0 0 0 1 ];

Good, that one is easy. The projection requires us to first calculate the focal length. Now, the field of view can be calculated by the width or height, and the focal length. Since the field of view is specified either by the VRML file, or defaulting to a standardised value, and the width or height is known, we can calculate the focal length. I won't go into calculating deriving this equation, because I don't know anything about trig identities. Thanks Wolfram|Alpha. Here is what you need:

v = 0.785398;
fx = 0.5 * rendW * cot(0.5 * v);
fy = 0.5 * rendH * cot(0.5 * v);

With this we can calculate the projection matrix.

K = [ fx 0  rendW/2 0
      0  fy rendH/2 0
      0  0  1       0
      0  0  0       1 ];

We need to add a column of ones to the points to make them homogeneous. This can be discarded later. Then, we apply these transformations to the points, and scale the X,Y values by Z.

P = K * C * [ P ones(size(P,1),1) ]';
P(1,:) = rendW - (P(1,:) ./ P(3,:));
P(2,:) = P(2,:) ./ P(3,:);
P = P';

That's it - you have the 2D projection your mesh. Unfortunately, during the projection, the Z axis was distorted. However, if you want to render these points as a cloud, you can, but you need to add a shift.

P(:,3) = P(:,3) + 500;
figure
imshow(I), hold on
plot3(P(:,1), P(:,2), P(:,3), '.', 'MarkerSize', 5)

And that's that.

Tags: math matlab hacks