How can you calculate the change in angle(RotX,RotY,Rotz) of a normal vector (say of that of a flat mirror) to a normal vector of an average plane thru a deformed mirror?

I need help with an issue concerning calculating the rotation angles given the normal vector of a mirror and the deformed normal vector of a mirror. I am able to find the rotation about the x-axis and y-axis but am having issues with calculating the rotation about the z-axis. I've attached the matlab script I wrote and have more details about the issue below.

NormAngleCalc.m can be summarized to the following steps/functions:

1. Finite element data points are inputted. The input data is comprised of the mirror [x,y,z] location points and the directional deformation [DeltaX, DeltaY, DeltaZ].

2. Normal vectors are calculated for the original mirror face and the deformed mirror face. For reference, the undeformed normal vector is parallel to the z-axis and the deformed normal vector is very close to being parallel to the z-axis.

3. A plane is calculated in the z-direction using the normal vectors from (2) and the center point of the mirror.

4. Three random points are determined on the plane which are used to calculate the normal of the plane.

5. The rotation angles to go from the undeformed to deformed normal vector are calculated.

Problem:

My goal is to use the original normal vector and the deformed normal vector to determine the rotation angles in the x-axis, y-axis, z-axis. Currently the script is working for rotations about the x-axis, y-axis, and rotations in the x-axis and y-axis simultaneously. For example, if I input Ansys data that was only rotated 2.5 degrees about the x-axis in the model, NormAngleCalc.m outputs the angle 2.5 degrees with the rotation axis [1 0 0].

The issue arises when I am trying to determine a rotation about the z-axis. For example, if I input Ansys data that was only rotated 2.5 degrees about the z-axis, the resulting output is an angle .0343 degrees with the rotation axis [.2174 -.9761 0]. The output is a rotation about the x-axis and y-axis only. How can I take the rotation axis and angle in the x,y axis and rewrite that as a rotation about the z-axis (if that's possible)? Or are there any ways I change change my formulas to get 2.5 degree angle in the [0 0 1] axis?

Code:

processData;
%determining center point deformation and Data

data = [XLocationm YLocationm ZLocationm XDirectionalDeformationm YDirectionalDeformationm ZDirectionalDeformationm];
%cp_x = -7.4923e-6; %from ansys node 13141
cp_x = -6.0007e-2; %from ansys node 13074
cp_y = 6.6728e-6; %from ansys node 13074
cp_z = -9.048e-3; %from ansys node 13074

xIndex= find(cp_x == data(:,1));
yIndex= find(cp_y == data(:,2));
zIndex= find(cp_z == data(:,3)); %z-loc is the same for the entire array

for i = 1: length(xIndex)
    for j = 1: length(yIndex)
        if xIndex(i) == yIndex(j)
            cpIndex = xIndex(i);
            break
        end
    end
end

cp_loc = [cp_x cp_y cp_z]; %[x y z] loc of center point
cp_def = [data(cpIndex, 4) data(cpIndex, 5) data(cpIndex,6)]; %deformation at [x y z]
cp_data = cp_loc + cp_def; %[x y z] loc + deformation
%% undeformed normal vector & plane calculation
V_ud = [XLocationm YLocationm ZLocationm];
V_ud = V_ud-mean(V_ud);
[u,s,w] = svd(V_ud,0);
norm_ud = w(end,:);

syms x_ud y_ud z_ud
P = [x_ud,y_ud,z_ud];
planefunc = dot(norm_ud,P-cp_loc);

z_plane_ud(x_ud,y_ud) = solve(planefunc,z_ud);
xLocMin = min(XLocationm);
xLocMax = max(XLocationm);
yLocMin = min(YLocationm);
yLocMax = max(YLocationm);

xRand_ud = (xLocMax-xLocMin)*rand(1,3)+xLocMin; %3 random x-points on plane
yRand_ud = (yLocMax-yLocMin)*rand(1,3)+yLocMin; %3 random y-points on plane
zRandOut_ud = double(z_plane_ud(xRand_ud,yRand_ud));
plane_xyz_ud = [xRand_ud' yRand_ud' zRandOut_ud'];
pAvg_ud = mean(plane_xyz_ud,1);
ckNorm_ud = null(plane_xyz_ud-pAvg_ud);

figure;
fimplicit3(planefunc); %y and z coords are switched --> needs to be fixed
hold on; 
quiver3(cp_x, cp_y,cp_z,norm_ud(1),norm_ud(2),norm_ud(3),'ShowArrowHead','off');
hold on;
quiver3(cp_x, cp_y,cp_z,ckNorm_ud(1),ckNorm_ud(2),ckNorm_ud(3),'ShowArrowHead','off');
legend('Plane','PlaneNorm','DataNormVector');
title('Undeformed Plane');
xlabel('x');
ylabel('y');
zlabel('z');

%% deformed normal vector & plane calculation

V_deformed = [xData yData zData];
V_deformed = V_deformed-mean(V_deformed);
[u_d,s_d,w_d] = svd(V_deformed,0);
norm_d = w_d(:,end);

syms x_d y_d z_d
P_d = [x_d,y_d,z_d];
planefunc_d = dot(norm_d,P_d-cp_data);

z_plane_d(x_d,y_d) = solve(planefunc_d,z_d);
xDataMin = min(xData);
xDataMax = max(xData);
yDataMin = min(yData);
yDataMax = max(yData);

xRand_d = (xDataMax-xDataMin)*rand(1,3)+xDataMin; %3 random x-points on plane
yRand_d = (yDataMax-yDataMin)*rand(1,3)+yDataMin; %3 random y-points on plane
zRandOut_d = double(z_plane_d(xRand_d,yRand_d));
plane_xyz_d = [xRand_d' yRand_d' zRandOut_d'];
pAvg_d = mean(plane_xyz_d,1);
ckNorm_d = null(plane_xyz_d-pAvg_d);

figure;
fimplicit3(planefunc_d);
hold on;
quiver3(cp_data(1), cp_data(2),cp_data(3),norm_d(1),norm_d(2),norm_d(3),'ShowArrowHead','off');
hold on;
quiver3(cp_x, cp_y,cp_z,ckNorm_ud(1),ckNorm_ud(2),ckNorm_ud(3),'ShowArrowHead','off');
legend('Plane','PlaneNorm','DataNormVector');
title('Deformed Plane');
xlabel('x');
ylabel('y');
zlabel('z');
%% plot deformed plane and normal vector

figure;
ezmesh(z_plane_d, [min(xData), max(xData), min(yData), max(yData), min(zData), max(zData)]);
hold on;
scatter3(xData,yData,zData,'MarkerEdgeColor','g')
hold on;
quiver3(cp_data(1),cp_data(2),cp_data(3),ckNorm_d(1),ckNorm_d(2),ckNorm_d(3),'ShowArrowHead','off');
% set(gca,'XLim',[-.1 .1],'YLim',[-.1 .1],'ZLim',[-.00904805 -.00904795]);
title('Deformed Plane and Normal Vector');
xlabel('x');
ylabel('y');
zlabel('z');
colorbar;
%% plot undeformed plane and normal vector
figure;
ezmesh(z_plane_ud, [min(XLocationm), max(XLocationm), min(YLocationm), max(YLocationm) min(ZLocationm) max(ZLocationm)])
hold on;
scatter3(XLocationm,YLocationm,ZLocationm,'MarkerEdgeColor','r')
hold on;
quiver3(cp_x, cp_y,cp_z,ckNorm_ud(1),ckNorm_ud(2),ckNorm_ud(3),'ShowArrowHead','off');
set(gca,'XLim',[-.1 .1],'YLim',[-.1 .1],'ZLim',[-.00905 -.009045]);
title('Undeformed Plane and Normal Vector');
xlabel('x');
ylabel('y');
zlabel('z');

%% determine angle difference between two vectors

norm_undeformed = ckNorm_ud;
norm_deformed = ckNorm_d;

%try 1
theta = rad2deg(atan2(norm(cross(norm_undeformed, norm_deformed)), dot(norm_undeformed,norm_deformed)));
theta_xyz = rad2deg(cross(norm_undeformed,norm_deformed));

%try 2
dir_cos_d = rad2deg(acos(norm_deformed./norm(norm_deformed)))';
dir_cos_ud = rad2deg(acos(norm_undeformed./norm(norm_undeformed)))';
ang_diff = dir_cos_d-dir_cos_ud;

%try 3
rotvec = vrrotvec(norm_undeformed,norm_deformed);
rotAxis = rotvec(1:3);
rotAng = rad2deg(rotvec(4));
rotAng_xyz = rotAng*rotAxis;
% 
% rotMat = vrrotvec2mat(rotvec);
% 
% Rx = rotx(rotAng_xyz(1));
% Ry = roty(rotAng_xyz(2));
% Rz = rotz(rotAng_xyz(3));

Answer 1

[Matt J]

The rotation angles to go from the undeformed to deformed normal vector are calculated.

There is no unique rotation that will take a single vector to a target vector. However, one particular choice is implemented by quadricFit.vecrot in this FEX download

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

R=quadricFit.vecrot(normalStart, normalTarget)

This gives you a 3x3 rotation matrix R for the rotation. There are 12 different ways to decompose this R into Euler angles (RotX,RotY,Rotz). You can use rotm2eul if you have the appropriate toolbox. If not, you can use the attached .m file rot2taitbryan whichgives a choice of 6 out of the 12 possible decompositions.