0

The title is a bit convoluted but I don't know a good way to express it. I am wondering whether this expression that I've worked out for my specific task has been discussed in the past since I have not found any mention of it, and most work tackling similar tasks don't go as far in their optimization.

My specific task is as follows: I have two vectors ReferenceBase and ReferenceTarget, which together implicitly describe a rotation that rotate the former onto the latter. I then have a third vector, just Vector, onto which I want to apply that same rotation. My first intuition was to compute the quaternion that represents this rotation and apply it with quaternion multiplication (using the quaternion(xyz = vector, w = scalar) convention):

def TransferUnitVectorRotation(ReferenceBase, ReferenceTarget, Vector):
    HalfVector = normalize(ReferenceBase + ReferenceTarget)
    Q = Quaternion(cross(ReferenceBase, HalfVector), dot(ReferenceBase, HalfVector)) # 9 mult
    return Q * Quaternion(Vector, 0) * Q.conjugate() # 32 mult

I then expanded and simplified the last line using properties of rotation quaternions to land on the well-known formula (that you can for example find here):

def TransferUnitVectorRotation(ReferenceBase, ReferenceTarget, Vector):
    HalfVector = normalize(ReferenceBase + ReferenceTarget)
    Q = Quaternion(cross(ReferenceBase, HalfVector), dot(ReferenceBase, HalfVector)) # 9 mult
    T = 2 * cross(Q.xyz, Vector) # 6 mult
    return Vector + Q.w * T + cross(Q.xyz, T) # 9 mult

But I then realized that I could expand and simplify even further and express everything in terms of ReferenceBase, HalfVector and Vector. Some more work involving the vector triple product led me to this expression that I have not found anywhere online:

def TransferUnitVectorRotation(ReferenceBase, ReferenceTarget, Vector):
    HalfVector = normalize(ReferenceBase + ReferenceTarget)
    HB = dot(HalfVector, ReferenceBase) # 3 mult
    VB = dot(Vector, ReferenceBase) # 3 mult
    VH = dot(Vector, HalfVector) # 3 mult
    return Vector + ((HB * VB * 2 - VH) * HalfVector - VB * ReferenceBase) * 2 # 7 mult

Bringing the total multiplication count from the original 41 to a mere 16 for the full operation (excluding the first line).

Now this looks an awful lot like an expression involving projections and symmetries, but I have had a hard time researching existing solutions to this task. The best I have found is topics such as this one, which only care about extracting a matrix or quaternion representing the rotation, when I am looking to then immediately apply said rotation to a given vector.

Here is some Python code testing the validity of the derivation.

import numpy as np

def normalize(v):
    return v / np.sqrt(np.dot(v, v))


def QuatConj(q):
    return np.append(-q[0:3], q[3])
def QuatMult(q1, q2):
    return np.append(q1[3] * q2[0:3] + q2[3] * q1[0:3] + np.cross(q1[0:3], q2[0:3]), q1[3] * q2[3] - np.dot(q1[0:3], q2[0:3]))


def TransferUnitVectorRotation(ReferenceBase, ReferenceTarget, Vector):
    HalfVector = normalize(ReferenceBase + ReferenceTarget)
    Rotor = np.append(np.cross(ReferenceBase, HalfVector), np.dot(ReferenceBase, HalfVector))
    return QuatMult(Rotor, QuatMult(np.append(Vector, 0), QuatConj(Rotor)))[0:3]


def TransferUnitVectorRotationFast(ReferenceBase, ReferenceTarget, Vector):
    HalfVector = normalize(ReferenceBase + ReferenceTarget)
    HB = np.dot(HalfVector, ReferenceBase)
    VB = np.dot(Vector, ReferenceBase)
    VH = np.dot(Vector, HalfVector)
    return Vector + 2 * ((HB * VB * 2 - VH) * HalfVector - VB * ReferenceBase)


b = normalize(np.random.rand(3))
t = normalize(np.random.rand(3))
v = np.random.rand(3) * 10
print(TransferUnitVectorRotation(b, t, v))
print(TransferUnitVectorRotationFast(b, t, v))

Has this been discussed anywhere? I didn't exactly do anything ground-breaking so I assume this formula surely is known, but I can't seem to find any mention of anything that looks like it.

Matrefeytontias
  • 782
  • 1
    Note that a single vector and its target are not sufficient to determine the rotation. If we impose the constraint that the axis of rotation is orthogonal to both the base and target vector, then the solution becomes unique; perhaps this is along the lines of what you had in mind – Ben Grossmann Dec 11 '23 at 19:04
  • Yes, in my case I am simply choosing the axis of rotation that happens to be orthogonal to both reference and target, as is shown in the construction of the quaternion. – Matrefeytontias Dec 11 '23 at 19:06
  • Well, whatever it is that you are choosing to do with this approach, you should be aware that this is a specific choice that you had made. If your goal is to simply find any rotation that matches base to target (as your description of the task seems to imply), you might be able to find a quicker approach that doesn't satisfy this additional constraint – Ben Grossmann Dec 11 '23 at 19:13
  • Another approach to the problem (which applies the same rotation as you have arrived at with your quaternion construction) is to use Rodrigues' rotation formula. In this case, you can take $$ \mathbf k = \frac{\mathbf v_{\text{Base}} \times \mathbf v_{\text{Target}}}{|\mathbf v_{\text{Base}} \times \mathbf v_{\text{Target}}|}, \quad \cos \theta = \frac{ \mathbf v_{\text{Base}} \cdot \mathbf v_{\text{Target}}} {|\mathbf v_{\text{Base}}|,|\mathbf v_{\text{Target}}|}, $$ and $\sin \theta = 1 - \cos^2\theta$. – Ben Grossmann Dec 11 '23 at 19:16
  • 1
    I suspect that this is ultimately equivalent to the formula you arrived at by simplifying a vector triple product – Ben Grossmann Dec 11 '23 at 19:19
  • I can definitely see it, this makes sense. Thanks for the insight. – Matrefeytontias Dec 11 '23 at 19:25
  • Out of curiosity, I'm wondering why you're only tracking multiplications. IIRC addition/subtraction is significantly different cost from multiplication on most architectures, but I'm kind of surprised to see addition operations entirely omitted. I do not keep up to speed on how modern architectures do it, either! Can you make a comment on this to help me, please? – rschwieb Dec 11 '23 at 19:37
  • T = 2 * cross(Q.xyz, Vector) # 6 mult seems like this would be 9 multiplications as well, accounting for 2* each entry and 6 for the cross product itself. – rschwieb Dec 11 '23 at 19:43
  • As far as I know, multiplication is several times slower than addition on most if not all hardware, and so in tasks like these multiplication is usually the metric we're interested in. In all honesty it still is a good idea to track addition, I just haven't done it. Regarding your second comment, multiplying by two is usually optimized by the compiler down to adding the value to itself and so I do not count it as a multiplication. – Matrefeytontias Dec 11 '23 at 19:45
  • 1
    I did not know about the "faster quaternion-vector multiplication," but it is amusing. However, in a lot of applications, one has to apply a single quaternion to many vectors rapidly. When this is the case, it makes more sense to pay a small fixed cost to convert the quaternion to a 3x3 matrix, and then you're only paying 9 multiplications per application to a vector after that. The lesson I learned was the superiority of quaternions really shines in composition, but matrices are still pretty efficient at applying the transformation. – rschwieb Dec 11 '23 at 19:54
  • From what i've been reading, multiplication 4x slower than addition sounds like a pretty common worst case. Apparently for some processors it can be as good as 2x or nearly 1x. – rschwieb Dec 11 '23 at 19:59
  • Applying one quaternion to a lot of vectors seems to be the common use case indeed. In my case, I'm deforming a mesh with the result of a physics simulation on reference tangent vectors, so I'm applying a different rotation to every vector. – Matrefeytontias Dec 12 '23 at 11:41
  • If the base/target $(b,v)$ are normalized vectors, then the rotation matrix (in the $bv$-plane) can be constructed as $$\eqalign{ \lambda &= \left({\tt1}+v^Tb\right)^{-1} \ M &= vb^T-bv^T \ R &= I + M + \lambda M^2 \ v &= Rb \ }$$ – greg Dec 12 '23 at 17:53
  • If $b,v$ are normalized vectors written as quaternions with zero real part, then the rotation quaternion can be constructed as $q=(b+v)b/|b+v|$ if you apply it as $q(-)q^{-1}$, and we aren't in an edge case like $b=-v$.) – rschwieb Dec 12 '23 at 20:29
  • I forgot to mention that the matrix formulation is valid for ${\mathbb R}^{n}$ not just for ${\mathbb R}^{3}$ – greg Dec 12 '23 at 22:31

0 Answers0