Quaternion to Euler angle converter?

I would like to write a conversion function from (unnormalized) quaternions to Euler angles representing a 3D rotation (from accelerometer data). Does anyone have such a routine already written and is willing to share it? I know I am being lazy but I am new to quaternions and my head is already spinning (pun intended).

Thank you! I was looking for the inverse of your function, to convert quaternions back to angles, but I see that it's not that difficult to write, and the explanations on your post were very helpful.

It is not clear how you want to actually use or calculate the Quat-->Euler transformation. There is a method which avoids all calculation; simply use Igor's built in Gizmo operations. If you setup a dummy Gizmo, then (as an example)

•ModifyGizmo setQuaternion={-0.234854,  0.326305,  0.00631929,  0.915585}
•GetGizmo curRotation;print GizmoEulerA,GizmoEulerB,GizmoEulerC
  32.6847  -36.4821  -11.8307

From this example, you should be able to write a simple wrapper function for the Quaternion input, that returns the global Euler angle outputs.

First note that I do not recommend using Euler angles.  There is a fair amount of literature that explains why using Euler angles is usually not a good idea.

However, if you need to convert from a quaternion to Euler angles and you have decided which formats you are using I'd get there in two steps:

1.  Convert your quaternion to a rotation matrix.  This is described e.g., here: https://automaticaddison.com/how-to-convert-a-quaternion-to-a-rotation-…

2.  Convert the rotation matrix to Euler angles.  This is described e.g., here:

http://eecs.qmul.ac.uk/~gslabaugh/publications/euler.pdf

Also, check out MatrixOP functions for quaternions: quatToEuler(qIn,mode) and quatToMatrix(qIn).

 

Thank you both to s.r.chinn and Igor.

s.r.chinn: I was not familiar with the capabilities of gizmo that you pointed out. Those are exactly the operations I was looking for. Thank you!

Igor: I think I am aware of the reasons not to use Euler angles in general. My purpose is to extract quaternion data into a form that I can use to create a plot. The links you provided will help me better understand quaternions and decide whether my intuition about how to visualize the rotations they represent is incorrect. Thanks!

I am getting back to this task after learning more about quaternions and gizmo. I am trying to implement the suggestion by s.r.chinn to use Gizmo's built-in rotation function,

•ModifyGizmo setQuaternion={-0.234854,  0.326305,  0.00631929,  0.915585}
•GetGizmo curRotation;print GizmoEulerA,GizmoEulerB,GizmoEulerC
  32.6847  -36.4821  -11.8307

I am unclear about (at least) two steps. Here is what I'd like to do.

I have a 3-column matrix time series, Accelerometer, of acceleration obtained from a wearable sensor. The columns are for Sensor(x, y, z), meaning the sensor's coordinate frame. There are about 200,000-500,000 rows (30-60 minutes of recording at 128 Hz). The sensor is strapped to the chest at some arbitrary angle.

I'd like to obtain a new Accelerometer_Rot time series that is a rotated version of Accelerometer, such that the columns of Accelerometer_Rot are aligned with x, y, z in external space. I'd like to apply arbitrary rotations around the Sensor-x, y, z axes and visualize the resulting Accelerometer_Rot wave so I can determine the desired rotation by visual inspection.

Based on s.r.chinn's suggestion, I generated the following code:

•newGizmo
•AppendToGizmo Accelerometer,name=path0
•ModifyGizmo displayLastObject
•AppendToGizmo axes=defaultAxes,name=axes0
•ModifyGizmo showAxisCue=1

The example commands that rotate the path, suggested in the Igor manual p II-428, are

ModifyGizmo setQuaternion={sin(pi/4),0,0,cos(pi/4)}   // Z up, X right, Y away
ModifyGizmo appendRotation={0,0,sin(pi/4),cos(pi/4)}

My questions start here:

1. what is the difference between the setQuaternion and the appendRotation commands?

2. I understood quaternions to contain {real, imaginary-x, imaginary-y, imaginary-z} components. Thus the appendrotation command above says, "Rotate around the y axis by pi/2 and around the z axis by pi/2". Is this correct?

3. Then what rotation is specified in the setQuaternion command above? What axis is intended by the first element specified?

4. After I have applied the rotation, I would like to recover the rotated 3-column wave. How do I access the resulting rotated wave? If I execute

•getGizmo objectNameList

I get S_ObjectNames, which contains: path0;axes0;

Is "path0" the rotated time series? How do I access it?

 

Update:

I got this to work. The gizmo approach was not a good match for what I wanted. Instead, I programmed all necessary calculations based on the clear tutorial I found here:

https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html

Here is the code in case it's of use to anyone else. The first function takes a point XYZ, a rotation angle, and a rotation axis, and creates a rotated version of that point. The second function takes a time series of XYZ points and calls the first function for each row, and generates a rotated version of the 3D time series.

 

//From https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html
Function rotate3Dpoint(pointVec1x3, thetaDeg, rotAxis)
    wave pointVec1x3    //1 row x 3 cols
    variable thetaDeg
    wave rotAxis

    variable thetaRad = thetaDeg * pi/180
    Variable r0, r1, r2, r3, s0, s1, s2, s3


    Make /O /N=4 pointQuat, rotQuat, rotQuatInv
    pointQuat = {0,pointVec1x3[0][0],pointVec1x3[0][1],pointVec1x3[0][2]}
    rotQuat[0] = cos(thetaRad/2)
    rotQuat[1] = rotAxis[0] * sin(thetaRad/2)
    rotQuat[2] = rotAxis[1] * sin(thetaRad/2)
    rotQuat[3] = rotAxis[2] * sin(thetaRad/2)

    rotQuatInv[0] = rotQuat[0]
    rotQuatInv[1] = -1 * rotQuat[1]
    rotQuatInv[2] = -1 * rotQuat[2]
    rotQuatInv[3] = -1 * rotQuat[3]

    Make /O /N=4 prod1, prod2
    r0 = pointQuat[0]
    r1 = pointQuat[1]
    r2 = pointQuat[2]
    r3 = pointQuat[3]

    s0 = rotQuat[0]
    s1 = rotQuat[1]
    s2 = rotQuat[2]
    s3 = rotQuat[3]

    prod1[0] = r0*s0 - r1*s1 - r2*s2 - r3*s3
    prod1[1] = r0*s1 + r1*s0 - r2*s3 + r3*s2
    prod1[2] = r0*s2 + r1*s3 + r2*s0 - r3*s1
    prod1[3] = r0*s3 - r1*s2 + r2*s1 + r3*s0

    r0 = rotQuatInv[0]
    r1 = rotQuatInv[1]
    r2 = rotQuatInv[2]
    r3 = rotQuatInv[3]

    s0 = prod1[0]
    s1 = prod1[1]
    s2 = prod1[2]
    s3 = prod1[3]

    prod2[0] = r0*s0 - r1*s1 - r2*s2 - r3*s3
    prod2[1] = r0*s1 + r1*s0 - r2*s3 + r3*s2
    prod2[2] = r0*s2 + r1*s3 + r2*s0 - r3*s1
    prod2[3] = r0*s3 - r1*s2 + r2*s1 + r3*s0

    Make /O /N=(1,3) W_pointRotated

    W_pointRotated[0][] = {{prod2[1]},{prod2[2]},{prod2[3]}}
    //W_pointRotated[0][0] = prod2[1]
    //W_pointRotated[0][1] = prod2[2]
    //W_pointRotated[0][2] = prod2[3]

End//rotate3Dpoint

//Rotates a 3D wave
//Calls rotate3Dpoint for each point in a 3D wave
//Inputs:
//1) time series of xyz 3D point coordinates as
//Nrowsx3col 2D wave, where N=number of points in (usually time)
//series. If you have 3 separate x, y, z waves, use concatenate
//2) Rotation angle in degrees.
//3) Rotation axis as a string
//set to x, y, or z for the desired rotation axis
//(Could change this function to accept an arbitrary angle but
//this need is so rare that I will keep the convenience of just
//specifying axis with a single letter)
//Output: waveXYZrot, a rotated version of the input wave
Function Rotate3Dwave(waveXYZ, thetaDeg, rotAxisStr)
    wave waveXYZ
    variable thetaDeg
    string rotAxisStr
   
    variable index, NPoints
    string wName, wNameRot
    wName = nameOfWave(waveXYZ)
    wNameRot = wName + "_Rot"
    duplicate /O waveXYZ, $wNameRot
    wave waveXYZrot = $wNameRot
    make /O /N=(1,3) pointXYZ
    make /O /N=3 rotAxis
    rotAxis[0] = stringMatch(rotAxisStr, "x")
    rotAxis[1] = stringMatch(rotAxisStr, "y")
    rotAxis[2] = stringMatch(rotAxisStr, "z")
   
    NPoints = DimSize(waveXYZ, 0)
    for (index = 0; index < NPoints; index += 1)
        pointXYZ = {{waveXYZ[index][0]},{waveXYZ[index][1]},{waveXYZ[index][2]}}
        rotate3Dpoint(pointXYZ, thetaDeg, rotAxis)
        wave W_pointRotated
        //waveXYZ[index][] = W_pointRotated
        waveXYZrot[index][0] = W_pointRotated[0][0]
        waveXYZrot[index][1] = W_pointRotated[0][1]
        waveXYZrot[index][2] = W_pointRotated[0][2]
    endFor

End //Rotate3Dwave