I’m spawning a trigger box relative to another object using local offset. The engine I’m using does not provide a native function to make that conversion (at least in a way accessible to me), so I am performing the conversion from local offset to world coordinates manually. The engine provides world rotation of objects in Euler angles, which is subject to gimbal lock, so I’m using the 3×3 rotation matrix from the object’s transform instead.

However, I still seem to be experiencing what I believe is gimbal lock even when using the 3×3 rotation matrix rather than Euler angles. Normally when I have an object only rotated on one or two axes, the trigger is spawned correctly (in the “push” position of the lever):

However, there are certain rotations where local to world space conversion gets confused and the trigger does not spawn in the right location.

The trigger spawns near the center of the lever rather than in the desired location. Those levers’ world rotations (from nearest to furthest from the camera) are:

1. X=0, Y=-36, Z=32 degrees
2. X=0, Y=-36, Z=-32 degrees
3. X=0, Y=36, Z=32 degrees

Even though I am not using Euler angles to compute the local->world transformation, it looks like this is still gimbal lock? Is that accurate, or is there another error that could cause this? The code I am using to make the transformation:

struct Quaternion {
    double w, x, y, z;
};

Quaternion matrixToQuaternion(const RE::NiMatrix3& m) {
    Quaternion q;
    double trace = m.entry[0][0] + m.entry[1][1] + m.entry[2][2];

    if (trace > 0) {
        double s = 0.5 / sqrt(trace + 1.0);
        q.w = 0.25 / s;
        q.x = (m.entry[2][1] - m.entry[1][2]) * s;
        q.y = (m.entry[0][2] - m.entry[2][0]) * s;
        q.z = (m.entry[1][0] - m.entry[0][1]) * s;
    } else {
        if (m.entry[0][0] > m.entry[1][1] && m.entry[0][0] > m.entry[2][2]) {
            double s = 2.0 * sqrt(1.0 + m.entry[0][0] - m.entry[1][1] - m.entry[2][2]);
            q.w = (m.entry[2][1] - m.entry[1][2]) / s;
            q.x = 0.25 * s;
            q.y = (m.entry[0][1] + m.entry[1][0]) / s;
            q.z = (m.entry[0][2] + m.entry[2][0]) / s;
        } else if (m.entry[1][1] > m.entry[2][2]) {
            double s = 2.0 * sqrt(1.0 + m.entry[1][1] - m.entry[0][0] - m.entry[2][2]);
            q.w = (m.entry[0][2] - m.entry[2][0]) / s;
            q.x = (m.entry[0][1] + m.entry[1][0]) / s;
            q.y = 0.25 * s;
            q.z = (m.entry[1][2] + m.entry[2][1]) / s;
        } else {
            double s = 2.0 * sqrt(1.0 + m.entry[2][2] - m.entry[0][0] - m.entry[1][1]);
            q.w = (m.entry[1][0] - m.entry[0][1]) / s;
            q.x = (m.entry[0][2] + m.entry[2][0]) / s;
            q.y = (m.entry[1][2] + m.entry[2][1]) / s;
            q.z = 0.25 * s;
        }
    }

    return q;
}

RE::NiPoint3 rotateQ(const Quaternion& q, const RE::NiPoint3& v) {
    RE::NiPoint3 u = {(float)q.x, (float)q.y, (float)q.z};
    RE::NiPoint3 uv = {u.y * v.z - u.z * v.y, u.z * v.x - u.x * v.z, u.x * v.y - u.y * v.x};
    RE::NiPoint3 uuv = {u.y * uv.z - u.z * uv.y, u.z * uv.x - u.x * uv.z, u.x * uv.y - u.y * uv.x};

    return v + ((uv * (float)q.w) + uuv) * 2.0f;
}

RE::NiPoint3 localToWorld(RE::NiPoint3 localOffset, RE::TESObjectREFR* relativeTo) {
    // Get object's rotation matrix and convert it to a quaternion for easier use
    RE::NiTransform transform;
    relativeTo->GetTransform(transform);
    Quaternion q = matrixToQuaternion(transform.rotate);

    RE::NiPoint3 rotatedOffset = rotateQ(q, localOffset);
    return relativeTo->GetPosition() + rotatedOffset;
}

For reference, I’m using the Creation Engine whose axial implementation (according to these docs) is:

1. Extrinsic, left-handed angles applied in the order of Z, Y, and then X
2. X+ is right, Y+ is forwards, Z+ is up

and I’ve noticed that even using an Euler angles transformation that applies the angles in the right order still hits this problem. Papyrus function adapted from here:

float[] function GetPosXYZRotateAroundRef(ObjectReference akOrigin, ObjectReference akObject, float fOffsetX, float fOffsetY, float fOffsetZ)
        ; flip angles since it's left-handed
        float fAngleX = -akOrigin.GetAngleX()
        float fAngleY = -akOrigin.GetAngleY()
        float fAngleZ = -akOrigin.GetAngleZ()

        float myOriginPosX = akOrigin.GetPositionX()
        float myOriginPosY = akOrigin.GetPositionY()
        float myOriginPosZ = akOrigin.GetPositionZ()

        float fInitialX = fOffsetX
        float fInitialY = fOffsetY
        float fInitialZ = fOffsetZ

        float fNewX
        float fNewY
        float fNewZ

        ;Objects in Skyrim are rotated in order of Z, Y, X, so we will do that here as well.

        ;Z-axis rotation matrix
        float fVectorX = fInitialX
        float fVectorY = fInitialY
        float fVectorZ = fInitialZ
        fNewX = (fVectorX * Math.cos(fAngleZ)) + (fVectorY * Math.sin(-fAngleZ)) + (fVectorZ * 0)
        fNewY = (fVectorX * Math.sin(fAngleZ)) + (fVectorY * Math.cos(fAngleZ)) + (fVectorZ * 0)
        fNewZ = (fVectorX * 0) + (fVectorY * 0) + (fVectorZ * 1)        

        ;Y-axis rotation matrix
        fVectorX = fNewX
        fVectorY = fNewY
        fVectorZ = fNewZ
        fNewX = (fVectorX * Math.cos(fAngleY)) + (fVectorY * 0) + (fVectorZ * Math.sin(fAngleY))
        fNewY = (fVectorX * 0) + (fVectorY * 1) + (fVectorZ * 0)
        fNewZ = (fVectorX * Math.sin(-fAngleY)) + (fVectorY * 0) + (fVectorZ * Math.cos(fAngleY))

        ;X-axis rotation matrix
        fVectorX = fNewX
        fVectorY = fNewY
        fVectorZ = fNewZ        
        fNewX = (fVectorX * 1) + (fVectorY * 0) + (fVectorZ * 0)
        fNewY = (fVectorX * 0) + (fVectorY * Math.cos(fAngleX)) + (fVectorZ * Math.sin(-fAngleX))
        fNewZ = (fVectorX * 0) + (fVectorY * Math.sin(fAngleX)) + (fVectorZ * Math.cos(fAngleX))

        ;Return result
        float[] fNewPos = new float[3]
        fNewPos[0] = fNewX + myOriginPosX
        fNewPos[1] = fNewY + myOriginPosY
        fNewPos[2] = fNewZ + myOriginPosZ
        return fNewPos
endFunction

I’m slightly skeptical that this is gimbal lock since it happens both with a rotation matrix transformation & Euler angles applied in the correct order, but I’m not an expert in linear algebra. Is this gimbal lock, or some other issue? Is there a way to work around it?