Averaging an array of rotations and using them to apply torque

I'm building a VR interaction system, and adding a throwing component.

It's set up so that I have a continually updating array of positions for the held object, and when I release the held object it calculates the delta positions for each of the time steps between those last 5 positions and averages them. This leaves me with a smooth throw vector for impulse force without any jitters introduced from an unintentional wrist flick on release. This system is implemented and works great.

I'm now trying to do the same thing with rotations and adding a torque force on release. However, with my current implementation it seems like it just picks a random axis to rotate on upon release.

I have my updating array of quaternions implemented the same way as the positions, no problem there, the issue is either in my quaternion averaging code, or my delta rotation calculation code:

The following is my method for calculating delta rotations. Note that rotations is the array of n rotations from the last frames, most recent at 0 and oldest at n

Quaternion[] deltaRotations = new Quaternion[rotations.Length - 1];

for (int i = 0; i < deltaRotations.Length; i++)
{
    deltaRotations[i] = Quaternion.Inverse(rotations[i]) * rotations[i+1];
}

The following is my method for averaging the delta rotations.

Quaternion avgRot = Quaternion.identity;

int quatCount = 0;

foreach (Quaternion quat in deltaRotations)
{
    quatCount ++;
    avgRot = Quaternion.Slerp(avgRot, quat, 1 / quatCount);
}

The resulting forces from Rigidbody.addTorque(avgRot, ForceMode.Impulse) set it rotating in a seemingly random directions.

Anything jumping out at you guys as blatantly wrong?

EDIT

Per Ruzihim's answer, I have converted the code to use angle / axis representations of rotations and average those. This nails the proper axis for rotation upon object's release, but seems to give me the same strength of rotation every time. Code for the new calculate rotation force method below, not that I'm passing the array of deltaRotations into this function rather than passing in the raw rotations array and calculating the delta rotations after the fact:

    Vector3 averageRotationsAngleAxis(Quaternion[] rotations) // returns an angle / axis representation of averaged rotations[]
    {
        Vector3[] deltaAxes = new Vector3[rotations.Length];

        for (int i = 0; i < rotations.Length; i++)
        {
            float angle;
            Vector3 axis;
            rotations[i].ToAngleAxis(out angle, out axis);

            deltaAxes[i] = angle * axis;
        }

        Vector3 returnVec = averageVector3(deltaAxes);

        return returnVec;
    }

Marking Ruzihim's answer as correct, when I Debug.Log the result of the throw rotations using his method they come out looking right, the problem must be somewhere in between that and the impulse torque

Solution

Using slerp to average quaternions can result in some unexpected behaviours. For instance, if you have 2 delta rotations, one of 180 degrees around the y axis and one of -180 degrees around the y axis, slerping from one to the other would never produce an identity rotation(a rotation of 0 degrees arund the y axis). it would actually always be a 180 degree rotation in some direction. Try it yourself:

Quaternion a = Quaternion.AngleAxis(180f, Vector3.up); 
Quaternion b = Quaternion.AngleAxis(-180f, Vector3.up); 
Quaternion c = Quaternion.Slerp(a, b, 0.5f)

float angle; 
Vector3 axis; 
c.ToAngleAxis(out angle, out axis); 
Debug.Log(angle); 
Debug.Log(axis);

So, instead, I would recommend averaging the axis & angle of rotation, specifically using the global axes of rotation. To get the delta rotations in global axes, start with global axis delta rotation * start rotation = next rotation and solve algebraically:

                globalAxesDelta * startRot = nextRot
globalAxesDelta * startRot * inv(startRot) = nextRot * inv(startRot)
                           globalAxesDelta = nextRot * inv(startRot)

So to calculate the rotation deltas using global axes, use deltaRotation = rotations[i+1] * Quaternion.Inverse(rotations[i]) instead of deltaRotations[i] = Quaternion.Inverse(rotations[i]) * rotations[i+1]. Then, convert them to axis/angle form using ToAngleAxis:

Vector3[] deltaAxes = new Vector3[rotations.Length - 1];

for (int i = 0; i < deltaRotations.Length; i++)
{
    Quaternion deltaRotation = rotations[i+1] * Quaternion.Inverse(rotations[i]);
    float angle;  
    Vector3 axis;
    deltaRotation.ToAngleAxis(out angle, out axis);
    deltaAxes[i] = angle * axis;
}

Then average these global axes:

Vector3 avgAxis = Vector3.zero;

foreach (Vector3 angleAxis in deltaAxes)
{
    avgAxis += (1f/deltaRotations.Length) * angleAxis;
}

Then you would apply the torque using something like AddTorque(avgAxis * someTorqueFactor) to apply it along global axes.