How can I calculate the pixel position in a figure-windows of a 3D object in perspective projection?

I have an animation of a complex scene of 3D objects (surfaces and patches) in an axis that has 'projection' set to 'perspective', while my cam(era)pos(ition) and cam(era)target change during the animation. I want to draw a trace of one of the object in the window, as a line connecting the center of this object in the subsequent animation frames. Therefore, I need to calculate the pixel position of that object, and draw that on a (invisible) static plane axis in the same figure window.

I am stuck in calculating the pixel position.

[az,el] = view(axis1);
fov = camva(axis1);
viewmatrix = viewmtx(az,el,fov,campos(axis1)+arbitrary_offset);
% arbitrary_offset is needed to get the correct projection
pixelposition = viewmatrix*[xyz;ones(1,size(xyz,2))];
% xyz is the position in the 3D frame
line(axis2,pixelposition(1)./pixelposition(4),pixelposition(2)./pixelposition(4));

In order to match the calculated pixelposition to the rendered scene, I have to modify the axis scaling of the 2nd axis, the axis limits, and an arbitrary_offset for the position used in the call to viewmtx. And still, it's not perfect, and the offset has to be changed when I resize the figure window.

What I would like is a function that takes the axis parameters of the 3D axis (position, aspect ratio, axis limits, camera parameters like position, target, upvector, viewangle) and the axis limits of the 2D overlay axis, and calculates the correct corresponding 4x4 transformation matrix.

How do I do that? Are there Matlab functions for this purpose available? Or can I request somehow the pixelposition of a point in a 3D axis directly through some hidden property?

[J-G van der Toorn]

I see a main obstacle in getting this solved in a robust way: displaying axes with either PlotBoxAspectRatio or DataAspectRatio set to manual, causes the axes to be resized within the figure, without updating axis' position property. In fact, not a single (observable) property of the scaled axis is changing when the containing figure is being resized. The Position property of the axis is the containing pox, not the actual one. This 3D-to-2D function can only work well when the actual axis position and dimension, as calculated under-water by the Matlab graphics engine, is available.

Answer 1

[Mike Garrity]

I'm afraid that it's rather more complex than just getting the view matrix.

Let's walk through an example in detail. I'll assume that you're in 3D, with a perspective projection, and 'axis vis3d'.

Create 2 axes with a line each. We'll make axis1 3D, and axis2 2D in pixels. Then we'll try to match the lines.

fpos = get(gcf,'Position');
axis1 = axes;
axis2 = axes('Position',[0 0 1 1], ...
             'Visible','off','HandleVisibility','off','PickableParts','none', ...
             'XLim',[1 fpos(3)],'YLim',[1 fpos(4)]);
l1 = line(nan,nan,nan,'LineWidth',5,'Color',[.75 .75 .75],'Parent',axis1);
l2 = line(nan,nan,nan,'LineWidth',1,'Color','red','Parent',axis2);

Give l1 random coordinates, and give axis1 a random view.

l1.XData = randn(1,5);
l1.YData = randn(1,5);
l1.ZData = randn(1,5);

view(axis1,360*rand,180*rand-90);
axis(axis1,'vis3d')
axis1.Projection = 'perspective';

Pull coordinates out of line 1 and append W.

xyzw = [l1.XData; l1.YData; l1.ZData; ones(size(l1.XData))];

We need 4 things to convert xyzw to pixels. - Model transform - View transform - Viewport - Projection transform

Compute model transform

xl = axis1.XLim;
yl = axis1.YLim;
zl = axis1.ZLim;
xscale = 1/diff(xl);
yscale = 1/diff(yl);
zscale = 1/diff(zl);
model_xfm = [xscale,      0,      0, -xl(1)*xscale; ...
                  0, yscale,      0, -yl(1)*yscale; ...
                  0,      0, zscale, -zl(1)*zscale; ...
                  0,      0,      0,             1];

Compute view transform

[az,el] = view(axis1);          
view_xfm = viewmtx(az,el);
view_offset = ((axis1.CameraPosition-axis1.CameraTarget).*[xscale,yscale,zscale] ...
            + [1/2 1/2 1/2]);
view_offset = view_xfm*[view_offset,1]';
view_xfm(1:3,4) = -view_offset(1:3);

Compute viewport and aspect ratio

old_units = axis1.Units;
axis1.Units = 'pixels';
viewport = axis1.Position;
axis2.Units = old_units;
ar = viewport(3)/viewport(4);

Compute projection transform

fov = axis1.CameraViewAngle;
tanfov = tand(fov/2);
n = .1;
f = 10;
r = tanfov * ar * n;
l = -r;
t = tanfov * n;
b = -t;

proj_xfm = [2*n/(r-l),         0,  (r+l)/(r-l),            0; ...
                    0, 2*n/(t-b),  (t+b)/(t-b),            0; ...
                    0,         0, -(f+n)/(f-n), -2*f*n/(f-n); ...
                    0,         0,           -1,           0];

Now w can compute pixel coordinates.

Multiply xyzw by the 3 transform matrices.

ndc = proj_xfm*view_xfm*model_xfm*xyzw;

% Go from normalized device coordinates [-1 to 1] to pixels.
dc = [.5*(1+ndc(1,:)./ndc(4,:)); .5*(1+ndc(2,:)./ndc(4,:))];
dc(1,:) = viewport(1) + viewport(3)*dc(1,:);
dc(2,:) = viewport(2) + viewport(4)*dc(2,:);

% Place the resulting coordinates on l2.
l2.XData = dc(1,:);
l2.YData = dc(2,:);
l2.ZData = zeros(1,size(dc,2));

As I said, I've made some simplifying assumptions here. In particular, if you're letting the PlotBoxAspectRatio and DataAspectRatio be autocalculated, then you need to account for the factors described on this page of the documentation. But hopefully this will be enough to get you unstuck.

[J-G van der Toorn]

Hi Mike,

Thank you for your extensive answer. It certainly gives insight in how to approach this properly.

With

 axis1 = axes('DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1]);

I get the axis ratios correct, and the transformations keeps working correctly. (I changed the line

 axis1.Units = old_units; % was: axis2.Units = old_units;

but that is not relevant.)

However, when I modify the CameraTarget of axis1, e.g.

 axis1.CameraTarget = [0 0 0];

then the pixel coordinate calculation fails, and the red line no longer overlaps the grey line, which surprises me, as the CameraTarget value is taken into account in the calculation of the view_offset. I suppose this has to do with the fact that the center of the axes is no longer in the center of the figure window. However, I did not manage to solve this by accounting for it in the calculation.

Furthermore, a change of the CameraUpVector

 axis1.CameraUpVector

also breaks the calculation. This seems harder to solve, but I think this should be done by modification of the model_xfm matrix, by multiplication by a rotation matrix that accounts for the rotation of the camera with respect to its view direction.

I am trying to get this done, but if you already have a solution at hand, I would be glad to hear it.

Best regards, Jan-Gerard

[Mike Garrity]

No, the CameraTarget and CameraUpVector go into the view matrix. My use of viewmtx is the issue there. It's a very old function, and it actually predates those properties. When it was introduced, there was only the azimuth & elevation.

[J-G van der Toorn]

I made some progress. I realized that the vector from CameraPosition to CameraTarget is already accounted for by the view parameters az and el. And that the offset should therefore not have the 1/2 terms. I changed some lines in your code, and now it works, at least for arbitrary CameraTarget and CameraPosition values. The lines

 view_offset = ((axis1.CameraPosition-axis1.CameraTarget).*[xscale,yscale,zscale] ...
             + [1/2 1/2 1/2]);
 view_offset = view_xfm*[view_offset,1]';

are replaced by

 view_offset = view_xfm*model_xfm*[axis1.CameraPosition 1]';

And I get correct results when I choose random CameraTarget and CameraPosition values. However, I get a mismatch when the CamUpVector differs from [0 0 1]. I'm working on that issue now.

[J-G van der Toorn]

Now I have something that almost works. I calculate the camera rotation by projecting the CameraUpVector of the original axis into the camera space. Then I calculate the angle that vector on the image plane, and apply that rotation to nodes of the the projected object in the image plane.

 % The CameraUpVector in the camera space:
 camupvec = proj_xfm*view_xfm*model_xfm*[axis1.CameraUpVector';0];

 % Its angle in the camera plane:
 rot = pi+atan2(camupvec(1)/camupvec(4),camupvec(2)/camupvec(4));
 % The corresponding rotation matrix:
 camrot = [cos(rot) -sin(rot);sin(rot) cos(rot)];

 % Go from normalized device coordinates [-1 to 1] to pixels, including the camera rotation:
 ndc1 = camrot*[ndc(1,:)./ndc(4,:);ndc(2,:)./ndc(4,:)];

 % The actual pixel coordinates of the object:
 dc = repmat(viewport(1:2)',1,size(ndc,2)) + diag(viewport(3:4))*(1+ndc1)/2;

This works pretty well. However, if you look well, you will still notice discrepancies between the object in the 3D axis and that in 2D.

Something is not perfect yet. However, I have the feeling that that has more to do with the 'underwater' calculation of the displayed 3D object, than with the code above.

Furthermore, I noticed that this code (the complete function) fails when the figure dimensions are far from square, which has probably to do with the DataAspectRatio and/or PlotBoxAspectRatio of the axes involved. That should be solvable.

It would be very nice if Matlab's graph3d toolbox included a function to convert to pixel coordinates, in a correct, robust, and easy manner. That would facilitate automated positioning of annotations to 3D figures, other than by using text(x,y,z,....);

[Mike Garrity]

Yes, we're hoping to introduce some new functions to help with this. The problem is that they will have to handle all of the funny edge cases. Some of the PlotBoxAspectRatio cases, in particular, can be rather hairy.

It looks like you're very close. My guess is that your adjustment needs to be applied before the proj_xfm, rather than after. The fourth row in the projection matrix changes angles differently at the front and the back of the box, and that looks like what you are seeing.

[J-G van der Toorn]

Hi Mike, I hade such a feeling too, and indeed, that did the trick. I will post my final version on Matlab Central as soon is it is complete and documented. The last part is now:

 % Calculate the camera rotation in the view-frame:
 camupvec = view_xfm*model_xfm*[axis1.CameraUpVector';0];
 rot = atan2(camupvec(1),camupvec(2));
 % Create a rotation matrix:
 camrot = blkdiag([cos(rot) -sin(rot);sin(rot) cos(rot)],1,1);

 % Now w can compute pixel coordinates.
 % Multiply xyzw by the 3 transform matrices.
 ndc = proj_xfm*camrot*view_xfm*model_xfm*xyzw;

 % Go from normalized device coordinates [-1 to 1] to pixels.
 dc = repmat(viewport(1:2)',1,size(ndc,2)) + ...
      diag(viewport(3:4))*(1+[ndc(1,:)./ndc(4,:);ndc(2,:)./ndc(4,:)])/2;

It looks perfect, except that it fails due to the hairy effects of the PlotBoxAspectRatio when the figure becomes taller than wide. Should be fixable, but not today. Thanks for your help.

[Mike Garrity]

That's great!

[Dmitry]

Dear J-G van der Toorn, I'm trying to solve very similar problem. I'll be glad if you post your final code on matlabcentral.

Answer 2

[Ulrich Reif]

ds2fig on FX does a perfect job.

[J-G van der Toorn]

Indeed. You can find it on https://nl.mathworks.com/matlabcentral/fileexchange/43896-3d-data-space-coordinates-to-normalized-figure-coordinates-conversion