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.
