Precision when rotating vectors

I have a bunch of polygons, each defined by the N vertices contained in a Nx3 array "points". Every polygon is on a different plane. I want to rotate them such that each polygon's vertices lie on a plane that is parallel to the X-Y plane. To do this, I'm using rotation matrices Rx and Ry to rotate the vectors defined by the polygon vertices by an angle calculated from the normal to each plane.

v2 = points(1,:)-points(2,:);
v1 = points(1,:)-points(3,:);
N = cross(v1,v2);

theta_x = -atan2(N(2),N(3));
theta_y = atan2(N(1),N(3));

Rx = [
        1               0               0
        0               cos(theta_x)    sin(theta_x)
        0              -sin(theta_x)    cos(theta_x)
     ];
Ry = [
        cos(theta_y)    0              -sin(theta_y)
        0               1               0
        sin(theta_y)    0               cos(theta_y)
     ];

points = points - repmat([0,0,points(1,3)],size(points,1),1);

points_xy = zeros(size(points));
points_xy(1,:) = points(1,:);
for nVec = 2:size(points,1)
    V = points(nVec,:)-points(1,:);
    V_rot = Rx*(Ry*V');
    points_xy(nVec,:) = V_rot;
end

I find that this code rotates the polygons so that they are nearly parallel to the x-y plane, but not as close to parallel as I'd like. I tested the code, generating 100 polygons, each on a random plane. After rotating with the code above, the z-coordinate of the vertices of these polygons ranged between 0 +/- 0.2, with a mean = 0.0021 and std. dev. = 0.0698. I'd like to see this z-coordinate to fall within something like z = 0 +/- 0.01 (vertices are separated by ~1-2 units). Why such a large range of z values here? Does this have something to do with the precision of matlab's trig functions?

Answer 1

[Roger Stafford]

I have not actually tried out your code, but based on my reading of it, it appears to have a serious flaw. For the sake of discussion suppose that your polygon is a triangle whose plane is already parallel with the x and y axes. That will give a zero for the two theta's. Then in the line

 "points = points - repmat([0,0,points(1,3)],size(points,1),1);"

the three vertices are all projected vertically so they lie within the x,y plane itself. Finally in the for-loop Rx and Ry are identity matrices and the triangle is simply moved without rotation in a direction so that vertex #1 would move to the origin, except that only vertices #2 and #3 are actually moved. Vertex #1 has remained behind where it was, and the triangle has now assumed a completely different shape. I am sure this is not what you intended.

I will make two suggestions. First, there is no need to do two rotations. All you need is to rotate your polygon about a line which is orthogonal to both the polygon's normal and the z-axis, which is to say, about their cross product, and the amount of rotation should be such that the polygon's normal rotates to become parallel to the z-axis.

The second suggestion has to do with determining the normal direction of a polygon. Your method depends entirely on the first three vertices. However, if you have a polygon whose vertices deviate to some degree from being coplanar, it would be best to use matlab's 'svd' (singular value decomposition) function to determine the best normal direction. First subtract the mean values of the x, y, and z coordinates for the vertices from each vertex and then apply the 'svd':

 [U,S,V] = svd(bsxfun(@minus,points,mean(points,1)),'econ');

Your "best" normal vector will be the unit vector V(:,3) which has the smallest singular value.

Answer 2

[Brian Hannan]

Hi Roger. Thanks for the reply. Actually, I did intend for the vectors corresponding to each of the vertices to be rotated one at a time while leaving vertex #1 in place. Defining my vectors from the arbitrarily-chosen vertex #1 to all other verts and rotating them one-by-one should get all points onto the same plane, which is parallel to the XY axis and offset by the z coordinate of vertex #1, V1z. Subtracting each point by [0,0,V1z] should then drop each polygon onto the x-y plane.

While it is true that this code will not work for polygons whose vertices may define a plane, in this case I defined them to be coplanar. What I'm puzzled by is why I can rotate random polygons with this code and get them "close" to parallel with the x-y plane, but not reliably better than ~8 deg. parallel it. If these rotation matrices and angles are in fact defined correctly, then I believe that this method, while less elegant, should get the same result that could be achieved by a single rotation.