transfer function of a geometry to match another identical geometry

I have 2 identical 3D spatial geometry (in a CAD software). I want to find a transfer function (sequences of transforms and rotations around x, y and z axes) to move one to merge and match the other one. Is there any code or program that takes the coordinate of several points from both source and destination geometries and calculate the required moves and rotations?

Answer 1

[David Goodmanson]

Hi Hamid,

Do you need the actual sequence of rotations around the x,y, and z axes, or do you just need a 3x3 matrix that gets the job done regardless of the actual sequence? If it is the latter, then here is a method.

Assume you have a set of N points p_full in the form of an Nx3 matrix, where the x,y,z coordinates read across and each row represents one point. Assume the same for another set of M points in q_full. p_full and q_full must have at least three points in common or you can't get anywhere. So let p and q each be nx3 matrices for n points in common between p_full and q_full, with the points in the same order for p and q. n has to be at least 3, and more is better if there is some slight inaccuracy in the coordinates. The following code transforms q to qnew, which falls on top of p.

% make up some data
p = rand(12,3) + [1 4 3];
S = [orth([2 1 5]'),null([2 1 5])];
q = (p-mean(p))*S + [3 0 6];
% -----------------------

% find the rotation and translation from that data 
pcm = mean(p);
qcm = mean(q);
R = (q-qcm)\(p-pcm);
R*R'                     % check that R is an orthogonal transformation, 3x3 identity
qnew = (q-qcm)*R + pcm;  % qnew should fall on top of points p

figure(1)
scatter3(p(:,1),p(:,2),p(:,3),'k')
hold on
scatter3(q(:,1),q(:,2),q(:,3),'r')
hold on
scatter3(qnew(:,1),qnew(:,2),qnew(:,3),'b')
hold off

The new blue qnew points fall on top of the old black p points so you can't see them. If you have extra points in q_full that don't correspond to anything in p, the same transformation

qnew_full = (q_full-qcm)*R + pcm;

gives all the transformed points.

[Hamid]

Thanks David, it was a nice code (although I needed to adjust it a little bit because + sign does add the vectors to each row of the matrices). However, I need to have the rotation angles around x, y and z, because I need to perform the transformation inside my CAD software which rotates around the axes one by one. The other issue is that the corresponding points in the source and destination are not precisely matched due to inevitable tolerance error of the software. Let’s say I need to find the best match (not perfect) by increasing the number of points from both source and destinations to make the translation more accurate. In a perfect situation, only three points would be enough for a 3D “registration”. Do you have any suggestion?

[David Goodmanson]

Yes, it's good to use as many matching points as possible (assuming no points are seriously off for some reason). Matlab solves for R by minimizing the least squares error which is probably the best you can do.

I'm not sure what your first sentence means exactly. Adding the vector to each row of the matrix was what was supposed to happen.

Does the CAD software use fixed-in-space axes and rotate the model (active transformation)? Can you stretch or shrink the model along the individual axes? R can be decomposed reasonably accurately into a product such as Rz*Rx*Rz, and basically exactly if the size change along the axes is allowed.

Answer 2

[Hamid]

Thanks for the comments. In the earlier versions of Matlab, I could not add matrices with different orders (e.g., p = rand(12,3) + [1 4 3];). As for the transfer function, I decided to go with pcregrigid command which does not require the correct order of points and it can optimize the transfer function corresponding to the point clouds of fixed and moving objects. However, this command is in the vision toolbox, so your code is the best alternative if someone does not access to that toolbox.