Euler 3D rotation between two vectors

Hello everyone,

I would like to obtain the Euler angles needed to rotate a vector u = (0,0,1) to a vector v, defined between an arbitrary point (x,y,z) and the origin (0,0,0). So v = (x,y,z). I am currently working with a right-handed 3d coordinate system, where X is pointing to the right side of the current reader, Y is pointing upwards and therefore Z is pointing outside the screen of the computer.

I calculate the angles needed to rotate by taking as a reference axis Y -> X -> Z. I have tried to use Euler angles directly, but somehow the objects that I rotate end up in a different location that they should be. Searching on the internet, I have found loads of information about "rotation matrixes" and "quartenions", but I am lost about how to apply them to my case. Any help would be appreciated :)

EDIT: After some tinkering and searching on internet I managed to obtain an answer for a XYZ coordinated system. Still not useful for me, but maybe other people can profit from it.

a = [0 0 1].';
b = [0 5 0].';
% Method described in
% https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
a = a./norm(a);
b = b./norm(b);
v = cross(a,b);
vx = [0 -v(3) v(2) ; v(3) 0 -v(1); -v(2) v(1) 0 ];
c = dot(a,b);
I = eye(3);
R=I+vx+vx^2*(1/(1+c));
R = round(R,5);

% From Rot matrix to euler coordinates, follows XYZ, described in:
% http://www.gregslabaugh.net/publications/euler.pdf
if (R(3,1) ~=1) && (R(3,1) ~=-1)
    theta_1 = -asin(R(3,1));
    theta_2 = pi-theta_1;
    chi_1 = atan2((R(3,2)/cos(theta_1)),(R(3,3)/cos(theta_1)));
    chi_2 = atan2((R(3,2)/cos(theta_2)),(R(3,3)/cos(theta_2)));
    phi_1 = atan2((R(2,1)/cos(theta_1)),(R(1,1)/cos(theta_1)));
    phi_2 = atan2((R(2,1)/cos(theta_2)),(R(1,1)/cos(theta_2)));
    theta = min(theta_1,theta_2);
    chi = min(chi_1,chi_2);
    phi = min(phi_1,phi_2);
else
    phi = 0;
    if R(3,1) == -1
        theta = pi/2;
        chi = phi+atan2(R(1,2),R(1,3));
    else
        theta = -pi/2;
        chi = -phi+atan2(-R(1,2),-R(1,3));
    end
end

theta = rad2deg(theta)
chi = rad2deg(chi)
phi = rad2deg(phi)

[darova]

Do you have a picture or drawing what you are trying to implement?

[Jesus Sanchez]

Yes, sorry for the lack of information. I want to create several dipoles tangential to a sphere surface. This is an example, taking into account only XZ and YZ planes, which I was a ble to do correctly as there is always one coordinate that its 0:

To do this I:

1. Create the geometrical form at (0,0,1).
2. Calculate the vector rotation from there to (x,y,z), which is where I want it to be. I do this by projecting the vector (x,y,z) in YZ and XZ planes, and calculating their angle regarding (0,0,1).
3. First, I calculate the rotation taking X as the rotating axis -> alfa is the angle between (0,0,1) and the unit vector (0,y,z). The I do the same regarding the Y axis.

As you can see in the figure above, this approach works when one or more coordinates are 0. However, when I try to use it in a different way:

The matlab figure is the positions that I am trying to achieve on the sphere. However, you can see that the rotations are being calculated too large at the 4 points at the corners. Looking reasons about this I stumbled up on rotation matrixes and quartenions, but I am trying to apply them and I don´t really know how. I dont want to use robotics toolbox, as future students that use my code will probably not have a license to use it.

Edit 1: please note that in the image, the z coordinate is not 0, altought it has the same value for evey point, as it is a cut on the XY plane. I have tried to create different cuts in other planes too, but the result is usually wrong.

Did I explain the problem well? Do you happen to know why this rotation is being wrongly made?

Edit 2: The code that I use to achieve the rotations in both X and Y axis as stated before is as follows. The + or - signgs used on the angles are to cope with the rotations within the right-hand coordinate system.

% 4. Rotate
    A = insec(i,:) - [0 0 0]; % Vector between spherical point and origin
    normA = norm(A);
    uA = A./normA; % Unity vector
    z = [0 0 1];
    uAx = [uA(1) 0 uA(3)];
    uAy = [0 uA(2) uA(3)];
    %- Rotation must be done axis per axis.
    %-- ORDER OF ROTATIONS: YX -> Z not used (yet). 
    % Right hand rotation system
%[...] Lots of code that handles special cases where coordinates are 0
...
    
else % Finally, no coordinate is 0. Octant-situational
           if uA(3) > 0 % Left side of 3D coordinate plane.
                if uA(2) > 0 % Positive Y coordinates
                    alpha(i) = -rad2deg(acos(z*uAy.'/(norm(z))*norm(uAy)));
                    if uA(1) > 0
                        beta(i) = rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    else
                        beta(i) = -rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    end
                else % Negative Y coordinates
                    alpha(i) = rad2deg(acos(z*uAy.'/(norm(z))*norm(uAy)));
                    if uA(1) > 0 % Does not change in this side
                        beta(i) = rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    else
                        beta(i) = -rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    end
                end
            elseif uA(3) < 0 % Right side of 3D coordinate plane
                if uA(2) > 0 % Positive Y coordinates
                    alpha(i) = -rad2deg(acos(z*uAy.'/(norm(z))*norm(uAy)));
                    if uA(1) > 0
                        beta(i) = -rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    else
                        beta(i) = rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    end
                else % Negative Y coordinates
                    alpha(i) = rad2deg(acos(z*uAy.'/(norm(z))*norm(uAy)));
                    z = [0 0 -1];
                    if uA(1) > 0
                        beta(i) = -rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)))-180;
                    else
                        beta(i) = rad2deg(acos(z*uAx.'/(norm(z))*norm(uAx)));
                    end
                end
            else
                error('Something wrong at octants');
            end
        end
    

Answer 1

[darova]

Here is rotation using Rodriguez matrix

Here is rotation about Y axis and Z.

Here is the difference between them

See attached script

[Jesus Sanchez]

But how can I transform the matrix in the angle that I need to input to my 3d program? I guess it has something to do with a 2-step axis-angle representation, using first the X axis and then the Y axis. But how can the rotation matrix be decomposed in that?

[darova]

I used dot product to calculate angle between vectors. I rotated the model about Oy axis

th = acos(dot(a,b));        % angle between vectors

Use atan2 to calculate the angle rotation about Oz axis

ph = atan2(b(2),b(1));      % angle to rotate about Oz axis

[Jesus Sanchez]

It works! I have to wait until monday to try it with the 3D simulation program, but at least the demo with matlab looks great. Thank you very much!

[Chad]

Dear All,

Thank you for putting this example together. I am close to fully understanding this example but have some questions. say I start with a vector at:

theta1=0
phi1=0

I get the x,y,z from this theta_1,phi_1

s1 = cosd(theta_1)*sind(phi_1)
s2 = sind(theta_1)*sind(phi_1)
s3 = cosd(phi_1)

Now I want to rotate to a new theta_2,phi_2

theta2=90
phi2=135
t1 = cosd(theta_2)*sind(phi_2)
t2 = sind(theta_2)*sind(phi_2)
t3 = cosd(phi_2)

a = [s1 s2 s3];
b = [t1 t2 t3];
a = a./norm(a);
b = b./norm(b);

Now the process to generate the Rotation Matrix

th = acos(dot(a,b));        % angle between vectors
u = cross(a,b);             % vector of rotation
u = u/norm(u);
K = [0 -u(3) u(2)
    u(3) 0 -u(1)
    -u(2) u(1) 0];
R = @(t) eye(3) + sin(t)*K + (1-cos(t))*K^2;    % rotation matrix

Here is my question. How do I obtain the "data to rotate" below?

x = [-1 1 0 0];         % data to rotate
y = [0 0 -1 1];
z = x*0 + 1;

Is there a way to generate this from s1,s2,s3 above?

Now the x,y,z is rotated.

v = R(th)*[x;y;z];          % rotate data
[x1,y1,z1] = deal( v(1,:),v(2,:),v(3,:) );

My second question is how do I convert x1,y1,z1 to my theta_2,phi_2?

Any help is appreaciated.

CJ