Rotation (preferable a quaternion) between two vectors

[Elroy]

This is probably mostly directed at The Trick, but any/all are certainly welcome.

I need some code that will return a rotation (quaternion, or transformation matrix if we must) between two vectors. Ideally, the vectors would be declared with D3DVECTOR, but you could use something as simple as the following (for the vectors and quaternion).

Public Type VecType
    ' This is used for both 3D spatial vectors, and Euler angles as well.
    X As Double
    Y As Double
    Z As Double
End Type

Public Type QuatType                      ' To get axis and theta back out (JPL convention):
    w As Double         ' cos(theta/2)              theta = ACos(q.w) * 2
    X As Double         ' v.x * sin(theta/2)        v.x = q.x / sin(theta/2)    solve for theta first, and then plug in here.
    Y As Double         ' v.y * sin(theta/2)        v.y = q.y / sin(theta/2)
    Z As Double         ' v.z * sin(theta/2)        v.z = q.z / sin(theta/2)
End Type

I've got some C code that does it I wrote a while back, and I'll probably be working on translating that code. I was just wondering if there was already something out there in VB6.

Just to say it again, imagine two 3D vectors emanating from some origin. There is some rotation that could be performed on one, such that it would be rotated parallel to the other vector. (It would actually perfectly overlay the other.) It's that rotation I'm after.

Public Function QuatBetweenVectors(V1 As D3DVECTOR, V2 As D3DVECTOR) As D3DQUATERNION
    QuatBetweenVectors = ????
End Function

--- or ---

Public Function QuatBetweenVectors(V1 As VecType, V2 As VecType) As QuatType
    QuatBetweenVectors = ????
End Function

Thanks In Advance,
Elroy

[Elroy]

Ok, here's my first pass at a translation. I've still got to test though, and I'm certainly open to critiques:

I just did it with UDTs so people wouldn't have to mess with any DirectX stuff.

Option Explicit
'
Public Type VecType
    ' This is used for both 3D spatial vectors, and Euler angles as well.
    X As Single
    Y As Single
    Z As Single
End Type
'
Public Type QuatType                      ' To get axis and theta back out (JPL convention):
    w As Single         ' cos(theta/2)              theta = ACos(q.w) * 2
    X As Single         ' v.x * sin(theta/2)        v.x = q.x / sin(theta/2)    solve for theta first, and then plug in here.
    Y As Single         ' v.y * sin(theta/2)        v.y = q.y / sin(theta/2)
    Z As Single         ' v.z * sin(theta/2)        v.z = q.z / sin(theta/2)
End Type
'


Public Function QuatBetweenVectors(uV1 As VecType, uV2 As VecType) As QuatType
    ' Vector uV1 is the start, vector uV2 is the end.
    ' So, to make the two vectors parallel, we would do: uV_New = VecByQuatMultiply(uV1, QuatBetweenVectors(uV1, uV2))
    '
    Dim fAABB As Single
    Dim fAB As Single
    Dim uV3 As VecType
    Dim fCC As Single
    Dim fS As Single
    Dim fM As Single
    '
    fAABB = Sqr(VecDot(uV1, uV1) * VecDot(uV2, uV2))     ' Product of the lengths of the arguments.
    If fAABB = 0! Then
        QuatBetweenVectors = QuatIdentity
        Exit Function
    End If
    '
    fAB = VecDot(uV1, uV2) / fAABB                      ' Normalized dot product of the arguments (cosine).
    uV3.X = (uV1.Y * uV2.Z - uV1.Z * uV2.Y) / fAABB     ' Normalized crossproduct of the arguments.
    uV3.Y = (uV1.Z * uV2.X - uV1.X * uV2.Z) / fAABB
    uV3.Z = (uV1.X * uV2.Y - uV1.Y * uV2.X) / fAABB
    '
    fCC = VecDot(uV3, uV3)                              ' Squared length of the normalized crossproduct (sine).
    If fCC Then                                         ' Test if the arguments are not parallel or anti-parallel.
        If fAB > -0.707106781186548 Then                 ' Sqr(2)/2.
            fS = 1! + fAB                               ' Use the cosine to adjust the fS-element.
        Else
            fS = fCC / (1! + Sqr(1! - fCC))             ' Use the sine to adjust the fS-element.
        End If
        fM = Sqr(fCC + fS * fS)                         ' The magnitude of the quaternion.
        '
        QuatBetweenVectors.w = fS / fM                  ' Return the normalized quaternion.
        QuatBetweenVectors.X = uV3.X / fM
        QuatBetweenVectors.Y = uV3.Y / fM
        QuatBetweenVectors.Z = uV3.Z / fM
    Else                                                ' Parallel or anti-parallel.
        If fAB > 0! Then                                ' The arguments are parallel.
            QuatBetweenVectors = QuatIdentity
        Else                                            ' The arguments are anti-parallel.
            fM = Sqr(uV1.X * uV1.X + uV1.Y * uV1.Y)     ' The length of one argument projected on the XY-plane.
            If fM <> 0! Then                            ' Return uV1 rotation with the axis in the XY-plane.
                QuatBetweenVectors.w = 0!               ' 180 degrees.
                QuatBetweenVectors.X = uV1.Y / fM
                QuatBetweenVectors.Y = -uV1.X / fM
                QuatBetweenVectors.Z = 0!
            Else                                        ' The arguments are parallel to the Z-axis, rotate around the X-axis.
                QuatBetweenVectors.w = 0!               ' 180 degrees.
                QuatBetweenVectors.X = 1!
                QuatBetweenVectors.Y = 0!
                QuatBetweenVectors.Z = 0!
            End If
        End If
    End If
End Function

Public Function VecCross(uV1 As VecType, uV2 As VecType) As VecType
    ' Determines the cross-product of two 3D vectors.
    VecCross.X = uV1.Y * uV2.Z - uV1.Z * uV2.Y
    VecCross.Y = uV1.Z * uV2.X - uV1.X * uV2.Z
    VecCross.Z = uV1.X * uV2.Y - uV1.Y * uV2.X
End Function


Public Function VecDot(uV1 As VecType, uV2 As VecType) As Single
    ' Determines the dot product of two 3D vectors.
    VecDot = uV1.X * uV2.X + uV1.Y * uV2.Y + uV1.Z * uV2.Z
End Function

Public Function QuatIdentity() As QuatType
    QuatIdentity.w = 1!
End Function

Best Regards,
Elroy

EDIT1: Ohhh, I did it all as single precision, which is fine for my needs. It could certainly all be converted to double precision though.

EDIT2: Ok, testing now done, and I found a bug. That 0.707106781186548 is suppose to be -0.707106781186548. I also fixed it in the above code so people wouldn't copy-paste the bug.

[Elroy]

And attached to this one is my little test project, in case anyone's interested.

Also, in the BAS module, there are several other functions I cobbled together for my testing. But it was still the QuatBetweenVectors function that's new, and I wanted.

I'll leave the thread open for a day or two, in case there are any other comments, but I actually feel pretty good about it now.

Enjoy,
Elroy

[The trick]

Hi Elroy. To find the rotation (as a quaternion) you just need to find the rotation axis (cross-product) and angle (dot-product), then you can construct quaternion based on that data.

[Elroy]

Hi Trick,

Yeah, that's right. But it's not quite that easy. We've got to deal with both the parallel case (or at least find it), and the anti-parallel case. And it's a good idea to check for nearly parallel, and handle that a bit differently too (basically an epsilon check). And also, I'm trying to keep all my quaternions normalized (sqr(q.w^2 + q.x^2 + q.y^2 + q.z^2)=1).

Hopefully, I've got all those basis covered in the above procedure. At least I think I do. Actually, the epsilon check wasn't quite right in the above code. I fixed that as well (just setting it to 0.000001. Here's a version of it that uses the UDTs out of your TypeLib:

Public Function DxQuatBetweenVectors(uV1 As D3DVECTOR, uV2 As D3DVECTOR) As D3DQUATERNION
    ' Vector uV1 is the start, vector uV2 is the end.
    ' So, to make the two vectors parallel, we would do: uV_New = VecByQuatMultiply(uV1, QuatBetweenVectors(uV1, uV2))
    '
    Dim fAABB As Single
    Dim fAB As Single
    Dim uV3 As D3DVECTOR
    Dim fCC As Single
    Dim fS As Single
    Dim fM As Single
    '
    fAABB = Sqr(DxVecDot(uV1, uV1) * DxVecDot(uV2, uV2))    ' Product of the lengths of the arguments.
    If Abs(fAABB) < 0.000001! Then                          ' The arguments are too small, return zero rotation
        DxQuatBetweenVectors = DxQuatIdentity
        Exit Function
    End If
    '
    fAB = DxVecDot(uV1, uV2) / fAABB                        ' Normalized dot product of the arguments (cosine).
    uV3.X = (uV1.Y * uV2.Z - uV1.Z * uV2.Y) / fAABB         ' Normalized crossproduct of the arguments.
    uV3.Y = (uV1.Z * uV2.X - uV1.X * uV2.Z) / fAABB
    uV3.Z = (uV1.X * uV2.Y - uV1.Y * uV2.X) / fAABB
    fCC = DxVecDot(uV3, uV3)                                ' Squared length of the normalized crossproduct (sine).
    '
    If fCC <> 0! Then                                       ' Test if the arguments are not parallel or anti-parallel.
        If fAB > -0.707106781186548 Then                    ' -Sqr(2)/2.
            fS = 1! + fAB                                   ' Use the cosine to adjust the fS-element.
        Else
            fS = fCC / (1! + Sqr(1! - fCC))                 ' Use the sine to adjust the fS-element.
        End If
        fM = Sqr(fCC + fS * fS)                             ' The magnitude of the quaternion.
        '
        DxQuatBetweenVectors.w = fS / fM                    ' Return the normalized quaternion.
        DxQuatBetweenVectors.X = uV3.X / fM
        DxQuatBetweenVectors.Y = uV3.Y / fM
        DxQuatBetweenVectors.Z = uV3.Z / fM
    Else                                                    ' Parallel or anti-parallel.
        If fAB > 0! Then                                    ' The arguments are parallel.
            DxQuatBetweenVectors = DxQuatIdentity
        Else                                                ' The arguments are anti-parallel.
            fM = Sqr(uV1.X * uV1.X + uV1.Y * uV1.Y)         ' The length of one argument projected on the XY-plane.
            If fM <> 0! Then                                ' Return uV1 rotation with the axis in the XY-plane.
                DxQuatBetweenVectors.w = 0!                 ' 180 degrees.
                DxQuatBetweenVectors.X = uV1.Y / fM
                DxQuatBetweenVectors.Y = -uV1.X / fM
                DxQuatBetweenVectors.Z = 0!
            Else                                            ' The arguments are parallel to the Z-axis, rotate around the X-axis.
                DxQuatBetweenVectors.w = 0!                 ' 180 degrees.
                DxQuatBetweenVectors.X = 1!
                DxQuatBetweenVectors.Y = 0!
                DxQuatBetweenVectors.Z = 0!
            End If
        End If
    End If
End Function

You are always an inspiration to me.

Thanks,
Elroy

EDIT1: And here are the couple of support functions for the above to run (in addition to The Trick's DX9 TypeLib).

Public Function DxVecDot(uV1 As D3DVECTOR, uV2 As D3DVECTOR) As Single
    ' Determines the dot product of two 3D vectors.
    DxVecDot = uV1.X * uV2.X + uV1.Y * uV2.Y + uV1.Z * uV2.Z
End Function

Public Function DxQuatIdentity() As D3DQUATERNION
    DxQuatIdentity.W = 1!
End Function

[The trick]

anti-parallel case

Do you mean two vectors which have opposite direction?

[Elroy]

Yes.