6

So I'm not that well versed in linear algebra so I'm struggling with this.

I have a unit vectors v. I want to find two angles(angle 1, rotation around x-axis, and angle 2, rotation around z-axis) such that when I rotate v by them it aligns the vector v with the y-axis. From this question I have a function that can find the angle between arbitrary vectors and returns a rotation. But this function returns 3 angles. Essentially there is an infinite number of 3d rotation that aligns v with the y-axis so I want the two unique angles.

This the code I have now, it requires numpy and scipy:

import numpy as np
import random
from scipy.spatial.transform import Rotation as R

def rotation_from_unit_vectors(a, b):
    v = np.cross(a, b)
    c = np.dot(a, b)
    s = np.linalg.norm(v)
    kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
    rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s ** 2))
    return  R.from_matrix(rotation_matrix)


y_axis = np.asarray([0.0, 1.0, 0.0])

alpha = random.uniform(0, 10)
beta = random.uniform(0, 10)
gamma = random.uniform(0, 10)

v = np.asarray([alpha, beta, gamma])
v = v / np.linalg.norm(v)

r = rotation_from_unit_vectors(v, y_axis)
print(r.as_euler('xyz', degrees = True))
print(r.apply(v))
  • Even those restricted rotations are not unique (even beyond the normal angular symmetries). Do you care which one you get? – Davis Herring May 15 '21 at 09:41
  • @DavisHerring Hmm can you give an example of two that are not unique beyond angular symmetries? I can't really picture it... –  May 15 '21 at 17:38
  • @DavisHerring For context: I want to describe the direction(So no magnitude ) of vector. Where a neutral direction would be aligned with the Y-axis. I thought that 2 rotations could uniquely identify such a direction. But from your comment you are saying this is not the case. Then I'm wondering: What does describe 3D direction uniquely with (preferably)2 parameters. –  May 15 '21 at 17:53
  • Starting from <1,0,1>, you can rotate either by -90 degrees to <1,1,0> and then by 45 to <0,1,0> or by 90 to <1,-1,0> and then by 135. Two rotations can identify every possible **orientation** (or rotation of a *featureless* object like a point or symmetric rod through the origin), but not injectively, so that such an orientation does not determine a unique rotation pair even of this restricted type. That might not matter for your use case. – Davis Herring May 15 '21 at 22:21
  • @DavisHerring Hmmm... That's right. However this problem would not occur if you can only do rotations from -90 degrees to +90 degrees right? The reason I want to have the unique rotations is that I want to later compare different directions and check if they are the same. Anyway I will check if your answer works soon and accept it! –  May 17 '21 at 13:19
  • I don’t think +/-90 is enough to transform all possible vectors. Comparing different directions can be done with a dot product and a tolerance. – Davis Herring May 17 '21 at 14:07
  • @DavisHerring I don't want to overstep. But is there a way to "invert" this rotation. E.g once this rotation is calculated find the rotation that aligns the y unit vector with the given vector? –  May 19 '21 at 11:34
  • Do you mean other than rotating back around *z* and **then** around *x* so you can just negate the angles? – Davis Herring May 19 '21 at 13:13
  • @DavisHerring Yes I mean still using the euler angles convention of first rotating around x and then around z. –  May 20 '21 at 07:33
  • Well, a very similar approach of working out how to get the right *z* (which must be the job of the rotation about *x*) and then the right *x* (which doesn’t have even symmetry like *y* does) will surely work. – Davis Herring May 20 '21 at 07:45

2 Answers2

2

Hum, non-standard problem, required thinking a little.

Given v1 and v2 you want to rotate_z(rotate_x(v1, alpha), beta) to be on the same direction as v2.

We know that the aligned vector can be obtained by simply scaling scaling v2, this will gives x1,y3,z3 = v3 = v2 * |v1| / |v2|. Since rotation around z-axis, does not affect the z coordinate, we can determine alpha such that the z coordinate of rotate_x(v1, alpha) equals z3. After that we determine the angle beta to align place the X and Y coordinates properly

import numpy as np
def alignment_angles(v1, v2):
    x1,y1,z1 = v1 # initial vector
    x2,y2,z2 = v2 # reference vector

    # magnitude of the two vectors
    r1 = np.sqrt(x1**2 + y1**2 + z1**2)
    r2 = np.sqrt(x2**2 + y2**2 + z2**2)
    # this will be the result when aligning v1 to v2
    # it has the magnitude of v1 and the direction of v2
    x3,y3,z3 = x2*r1/r2, y2*r1/r2, z2*r1/r2
    # the rotation in x must set the z coordinate to the
    # final value, since the rotation over the z axis will
    # not affect the z coordinate (this have two solutions)
    rho1 = np.sqrt(z1**2 + y1**2)
    if(abs(z3 / rho1) > 1):
        raise ValueError('Cannot align these vectors')
    alpha = np.arcsin(z3 / rho1) - np.arctan2(z1, y1);

    # apply the rotation to make easier to calcualte the next stage
    y1, z1 = (y1 * np.cos(alpha) - z1 * np.sin(alpha), 
              y1 * np.sin(alpha) + z1 * np.cos(alpha))
    np.allclose(rho1, np.sqrt(z1**2 + y1**2))
    #assert(np.allclose(z1, z3))
    # now it is just a matter of aligning (x1, y1) to (x3, y3)
    beta = np.arctan2(y3, x3) - np.arctan2(y1, x1)
    x1, y1 = (x1 * np.cos(beta) - y1 * np.sin(beta), 
              x1 * np.sin(beta) + y1 * np.cos(beta))
    
    # ensure the fotated v1 was correctly aligned
    assert(np.allclose([x1, y1, z1], [x3, y3, z3]))
    
    return alpha, beta

Then you just call

alignment_angles((1,2,3), (3,4,5))

or you can also use numpy arrays with 3 rows.

Initially I thought it would be an application of spherical coordinates, that would be the case if the axis for the second rotation was the z-axis rotated accordingly to the first rotation.

Edit

There are some vectors that cannot be aligned with a rotation on x and a rotation on y. Suppose you want to align the vector v1 = (1, 0, 0) to the vector v2 = (0, 0, 1) the rotation in x will not affect v1, it will always point in the direction x, then when you rotate around the z axis it will always be on the XY plan.

The example you gave was really giving the wrong answer because asin is not injective.

I changed the function to raise a value error when you cannot align the given vectors.

Bob
  • 13,867
  • 1
  • 5
  • 27
  • Thanks for your answer. It seems to work in many cases. But also in many cases the assert fails: For example: `alignemnt_angles((-0.017868310075997075, -0.3998056402886916, -0.008330974849768172), (0,1,0))`. I printed the assert values and the `z` does not match. `[1.249000902703301e-16, 0.39994466064372136, -0.016658333534559422], [0.0, 0.4002914334001792, 0.0]` –  May 13 '21 at 10:55
2

Taking advantage of the fixed target alignment, this can be done in a straightforward manner with just trigonometry:

import math

def to_y(x,y,z):
  rx=-math.atan2(z,y)              # or +math.atan2(z,-y)
  y2=y*math.cos(rx)-z*math.sin(rx) # -> (x,y2,0)
  return rx,math.atan2(x,y2)

The rotations are defined as counterclockwise when looking at the origin from +x or +z (the right-hand rule); the rotation direction is always that with the smaller magnitude, but it may be possible to find a physically smaller rotation as indicated in the comment. Note that the input need not be normalized, and NaN is never produced (unless it appears in the input).

Davis Herring
  • 36,443
  • 4
  • 48
  • 76