Equations for Squad? - Githubissues

Thanks for providing this great module!

I have recently worked quite a bit with quaternion splines, and I was wondering about the implementation of Squad, especially the handling of non-uniform parameter intervals:

https://github.com/moble/quaternionic/blob/074b626d0c63aa78479ff04ed41638931ca6693a/quaternionic/interpolation.py#L174-L187

Is the derivation of this available somewhere?

Because I have derived equations for (non-uniform) Squad recently (https://splines.readthedocs.io/en/0.3.0/rotation/squad.html#Non-Uniform-Parameterization) and I arrived at slightly different equations.

I might of course have made a mistake, but since my derivations are all available, it should be easy to pinpoint the error if that's the case.

The problem is that if I am correct, the Squad implementation behaves worse than my alternative implementation using De Casteljau's algorithm and Catmull-Rom-like tangents. A comparison between the two can be seen in the animations on the page I mentioned above.

I detailed my derivation in Appendix A.4 of this paper. Unfortunately, I don't recall all the subtleties because it's been a decade since I really thought about this — I basically copied this package's code from this code via this code (and have since gone on to copy it over to this code). But I think the main points are that (1) we can require continuity of the derivative across quaternion "keys", which gives us equation A36, and (2) we can set that derivative to be equal to the average of the derivatives that you would have on either side if you were doing linear interpolation.

Now, (2) is just an arbitrary condition that gives you a tractable result. Though I think it's a reasonable condition, there is certainly room for argument. It might not be very natural, and may lead to spurious accelerations around the keys. It might make more sense to require that the second derivative be continuous. (It would surely lead to much more complicated expressions for $A_i$ and $B_i$.) I'm not sure about your condition; my intuition about those time-step factors at the beginnings of the exponentials suggest that they push the velocities in the wrong directions. But that's just a guess.

I'm not sure what you mean when you say that the "Squad implementation behaves worse...". How do you define "worse"? One possibly useful measure would be to plot the magnitude of the angular velocity of the interpolant as a function of time. I would think that a "good" interpolant would have a smooth plot. I'd be interested to see what this looks like for your version of squad, my version of squad, and your CatmullRom method.

Incidentally, for historical reasons, you may be interested in Shoemake (1985), which I believe is Shoemake's original introduction of something like squad — though he didn't call it that in this reference.

Okay, you piqued my curiosity, so I went ahead and did the angular-velocity comparison I mentioned above. I've come to two main conclusions:

The splines methods do not actually reproduce the input data. I think your code is essentially doing something that I call "unflipping" the quaternions which can be subtly wrong when dealing with continuous-in-time functions relevant for interpolation.
The result from splines.Squad is definitely not natural, splines.CatmullRom is reasonable, but I would argue that quaternionic.squad is slightly better. So my conclusion is that those factors of the time-step sizes in the exponentials in your expressions are causing problems.

I took the example from your documentation, tweaked the times to make things a little more extreme, evaluated each spline on a very fine grid, and then plotted the magnitude of the angular velocity. I actually had to use a log scale because the angular velocity for splines.Squad goes totally crazy — like 30 times larger than the angular velocity for the others. The times at which we have "key" quaternions are shown as vertical dashed lines.

angular_velocities

All the code and some more commentary are here. (I'm not sure why saving the animations doesn't work for me, but they evaluate just fine in a live notebook.)

Thanks for your nice responses!

I detailed my derivation in Appendix A.4 of this paper.

Thanks for the paper, that sounds interesting. My mathematical power isn't very strong, so I haven't fully understood all derivations, but at least some parts look familiar. I'll have another look when I have more time.

But I think the main points are that (1) we can require continuity of the derivative across quaternion "keys", which gives us equation A36, and (2) we can set that derivative to be equal to the average of the derivatives that you would have on either side if you were doing linear interpolation.

Now, (2) is just an arbitrary condition that gives you a tractable result. Though I think it's a reasonable condition, there is certainly room for argument. It might not be very natural, and may lead to spurious accelerations around the keys.

Yeah, I think (2) is where our approaches differ.

Of course there is no absolute "right" or "wrong", but your approach sounds like what I showed as MeanVelocity (but I'm not sure if it is exactly the same) and I called "wrong" here: https://splines.readthedocs.io/en/0.3.0/euclidean/catmull-rom-properties.html#Wrong-Tangent-Vectors Just to avoid misunderstandings: this shows the Euclidean case, but I assume that the behavior will (at least to some degree) also be relevant for rotations.

It might make more sense to require that the second derivative be continuous. (It would surely lead to much more complicated expressions for Ai and Bi.)

Maybe, but that's an entirely different topic. We would lose local control or interpolation when we do that.

I would like to have "natural" rotation splines in my collection just for completeness' sake, but I'm not really interested in splines without local control.

I'm not sure about your condition; my intuition about those time-step factors at the beginnings of the exponentials suggest that they push the velocities in the wrong directions. But that's just a guess.

I haven't done my derivation in the quaternion domain, but I have derived the factors for the Euclidean case and then just naively used them also for the rotation splines.

I'm using the approach from Catmull & Rom (1974) with Lagrange interpolation followed by B-Spline blending, which I have straightforwardly extended to the non-uniform case. Using the Barry-Goldman algorithm (https://splines.readthedocs.io/en/0.3.0/euclidean/catmull-rom-barry-goldman.html) yields the exact same result.

The equations for Bézier control points (which are the starting point for splines.quaternion.CatmullRom) are derived here: https://splines.readthedocs.io/en/0.3.0/euclidean/catmull-rom-non-uniform.html#Using-Non-Uniform-B%C3%A9zier-Segments

The equations for quadrangle points (which are the starting point for splines.quaternion.Squad) are derived here: https://splines.readthedocs.io/en/0.3.0/euclidean/catmull-rom-non-uniform.html#Using-Non-Uniform-Quadrangle-Interpolation

In the Euclidean case, both lead to identical interpolants.

I'm not sure what you mean when you say that the "Squad implementation behaves worse...". How do you define "worse"?

That's a good question! I meant the erratic behavior in the last animation there: https://splines.readthedocs.io/en/0.3.0/rotation/squad.html#Non-Uniform-Parameterization When changing the time values a bit (as you have done in your example), it gets even worse. You called it "crazy", so I guess we are on the same page?

One possibly useful measure would be to plot the magnitude of the angular velocity of the interpolant as a function of time. I would think that a "good" interpolant would have a smooth plot. I'd be interested to see what this looks like for your version of squad, my version of squad, and your CatmullRom method.

I have actually done that before reaching out (and you have just done it in your second comment).

I was initially surprised that your quaternionic.squad was closer to my splines.quaternion.CatmullRom and much "better" than my splines.quaternion.Squad even though I consider yours "wrong".

But in the meantime, I got an idea for an hypothesis, see below ...

Incidentally, for historical reasons, you may be interested in Shoemake (1985), which I believe is Shoemake's original introduction of something like squad — though he didn't call it that in this reference.

Yes, I'm aware of that paper.

It's a common misconception to associate it with Squad, it has nothing to do with it. BTW, this is also wrong in your documentation https://moble.github.io/Quaternionic.jl/dev/functions_of_time/#Quaternionic.squad-Union{Tuple{T},%20Tuple{AbstractArray{Rotor{T},%201},%20AbstractVector{%3C:Real},%20AbstractVector{%3C:Real}}}%20where%20T.

It describes how to use De Casteljau's algorithm with Slerp, which I have regurgitated here: https://splines.readthedocs.io/en/0.3.0/rotation/de-casteljau.html

It then suggests a way how to choose the control quaternions, which I have reproduced here: https://splines.readthedocs.io/en/0.3.0/rotation/catmull-rom-uniform.html#Shoemake's-Approach

This approach is only used for the uniform case. I don't know a way how to extend this to the non-uniform case.

In my splines.quaternion.CatmullRom I followed the first step of this paper (the part about De Casteljau's algorithm), but for the second step I used a different approach, which I very clumsily and un-mathematically described here: https://splines.readthedocs.io/en/0.3.0/rotation/catmull-rom-uniform.html

And contrary to Shoemake's second part, I could extend this to non-uniform intervals: https://splines.readthedocs.io/en/0.3.0/rotation/catmull-rom-non-uniform.html

Okay, you piqued my curiosity, so I went ahead and did the angular-velocity comparison I mentioned above.

That's great!

The splines methods do not actually reproduce the input data. I think your code is essentially doing something that I call "unflipping" the quaternions

Yes, that's what I call "canonicalization": https://splines.readthedocs.io/en/0.3.0/rotation/quaternions.html#Canonicalization

I don't allow rotations bigger than 180 degrees and I automatically call splines.quaternion.canonicalized() on construction.

Whether or not allowing that would be a good idea is a whole separate discussion, I wouldn't want to get into this right now.

BTW: you might also be interested in this discussion, which is somewhat related: https://github.com/scipy/scipy/pull/17334

The result from splines.Squad is definitely not natural,

Yes, and the question is whether that's an intrinsic limitation or a programming mistake of mine.

splines.CatmullRom is reasonable, but I would argue that quaternionic.squad is slightly better.

Well, that's what we still have to find out.

From your angular velocity plots and the animations I'm not yet convinced either way.

But I have an alternative argument which I'd like to bring forward:

In your equations (A38a) and (A38b) for the quadrangle points, isn't it strange that the contribution by one "chord" is weighted by some time intervals, while the contribution from the other chord is not weighted? To me, this seems very suspicious: changing the chord will change the quadrangle point, but changing the time interval won't?

So my conclusion is that those factors of the time-step sizes in the exponentials in your expressions are causing problems.

They might well cause problems, but they might also be "correct"?

I have a hypothesis: the factors in my Squad are "correct", but because of the extreme time values, the quadrangle points are very far away (maybe on the other side of the unit hypersphere, maybe even wrapping around?) and therefore lead to severe errors in the interpolated values. The Bézier control points of splines.quaternion.CatmullRom on the other hand are not that far away, leading to fewer errors.

I guess this hypothesis is easy to test by printing the angles between the quaternions of the control polygon, but this will have to wait another day, because this answer has taken me already a long time ...

But what do you think, does that hypothesis sound reasonable?

Or do you have another hypothesis?

All the code and some more commentary are here.

Very nice, thanks for sharing this!

(I'm not sure why saving the animations doesn't work for me, but they evaluate just fine in a live notebook.)

I don't know, but it seems to work on nbviewer: https://nbviewer.org/urls/gist.github.com/moble/5d0f7cdffc90e30c27c374f9ccf107b0/raw/041ced784313e2b28f92708721a17d0208d0338a/Compare.ipynb

A little update on my hypothesis: I have checked the angles between control points and they indeed seem to be much bigger for the quadrangle points. See the bottom of this page: https://splines.readthedocs.io/en/wip/rotation/squad.html

Here's a copy of the values, first for splines.quaternion.CatmullRom:

[[np.degrees(q1.rotation_to(q2).angle) for q1, q2 in zip(s, s[1:])]
 for s in cr3.segments]

[[17.022364436790564, 135.06479911384824, 49.77767008520975],
 [29.866602051125813, 90.45178400459615, 40.38373450835629],
 [71.7933057926333, 96.10442436530326, 37.85425435739818],
 [70.97672692012152, 134.5273935255376, 68.5666275352161],
 [22.85554251173871, 57.18209204731734, 11.348242957860345]]

And then for splines.quaternion.Squad:

[[np.degrees(q1.rotation_to(q2).angle) for q1, q2 in zip(s, s[1:])]
 for s in sq3.segments]

[[30.063924306160327, 155.95310872103568, 36.475607327942335],
 [101.32113146650647, 306.7249656389471, 80.19651207109983],
 [25.374677647496483, 279.10261645462316, 105.5506389482895],
 [30.023292856402406, 105.33337836375514, 11.549888167817787],
 [103.94899351035986, 256.9689694401947, 67.64382968886072]]

This of course proves nothing, but I think it supports my hypothesis. The angles between quadrangle points are up to more than 300 degrees, that's close to half a turn around the unit hypersphere. This example isn't extreme, I'm sure that for more extreme time values, the angles would even wrap around, which would maybe explain the erratic behavior?

And an update about my claim of quaternionic.squad being "wrong":

As I said above, there is no absolute "right" or "wrong".

What I'm interested in, are centripetal splines (https://splines.readthedocs.io/en/0.3.0/euclidean/catmull-rom-properties.html#Centripetal-Parameterization), because of their property of guaranteeing that there are no cusps and self-intersections (which I really hope also extends to rotation splines! But I think it does ... at least the self-intersection goes away in this example: https://splines.readthedocs.io/en/0.3.0/rotation/catmull-rom-non-uniform.html#Parameterization).

Another, related (or maybe just another aspect of the same) property is that the curve always moves away from the previous key frame and towards the next one. And this is something I can actually test!

So I have created this gist (based on yours): https://gist.github.com/mgeier/575063c29c0b53781e128b49c3855382 (nbviewer).

This shows -- regarding centripetal-ness -- no indication that splines.quaternion.CatmullRom is "wrong", but it shows that quaternionic.squad is "wrong".

Interestingly, it also shows that splines.quaternion.Squad is "extremely wrong"!

Now the question is: is there any way to change any of the "squad" implementations to fulfill the centripetal properties?

To answer my own question:

is there any way to change any of the "squad" implementations to fulfill the centripetal properties?

Yes! It turns out that I just have to use my factors in quaternionic.squad, and voilà, it ceases to be "wrong". You can try it out with #42.

This might mean that my hypothesis above is probably wrong, and something else has to be fixed in my splines.quaternion.Squad. I'll try that another day.

I couldn't help but investigate the problem with splines.quaternion.Squad:

And it turns out that I had simply swapped the equations for the incoming and outgoing quadrangle point!

I have fixed this in https://github.com/AudioSceneDescriptionFormat/splines/commit/d21fea56ab6bd323beaaf8d00f43bacb0131768f (in the wip branch) and now it also passes the "centripetal" test in https://gist.github.com/mgeier/575063c29c0b53781e128b49c3855382.

Now its behavior is a lot less crazy than before, but it still becomes erratic when using extreme time values. The hypothesis I mentioned above might still turn out to be true?

I have updated the text at https://splines.readthedocs.io/en/wip/rotation/squad.html#Non-Uniform-Parameterization with my new findings.

What's your opinion on all that?

From everything you've written, I wonder if one issue is that we just have very different motivations. It seems like you're more interested in the case where someone is coming in and designing motion by hand, manually selecting orientations, whereas I'd be more interested in physical motions that are actually being measured (in some broad sense of that word), and we have to approximate the motion between those measurements.

My reasons for thinking this are

You assume that the first and last quaternions are identical.
You

don't allow rotations bigger than 180 degrees

and are disturbed by quaternions that go to

the other side of the unit hypersphere, maybe even wrapping around?
You are concerned about "centripetal-ness", and want to avoid cusps and self-intersections.

None of these are important to me, and in particular that last one is actually undesirable to me. I generally deal with problems where the quaternions and the times at which they are measured are already known; we don't get to reparameterize.

It's also important to recognize that the two-dimensional analogy won't work very well in the three-dimensional space of unit quaternions. Cusps and self-intersection will be much more unusual, just for geometric reasons. But even when they do occur, for me that will be because some physical system actually behaved that way.

splines.CatmullRom is reasonable, but I would argue that quaternionic.squad is slightly better.

Well, that's what we still have to find out.

Obviously, it depends on your criteria, but I feel like CatmullRom starting off in the wrong direction is a bad sign, as is its larger discontinuity in the magnitude of the angular velocity. And your (corrected) version of Squad has this totally anomalous acceleration just after the

Another criterion is minization of "torque", which I measure as the average magnitude of the angular acceleration. For the example from my gist, CatmullRom's average is 27.86, your new version of Squad's average is 23.00, while quaternionic.squad's is 19.68. Now, that's just for one example. Interestingly, if I generate 5,000 totally random unit quaternions, and 5,000 random time intervals between 0.05 and 1, CatmullRom actually has somewhat lower torque on average compared to quaternionic.squad, while your new version of Squad still has significantly more. I wonder if the bad results for the simple case were just because of problems with the treatment of endpoints.

I've updated my gist with your suggested revision replacing the buggy Squad, and showing these torques. My conclusion is still that quaternionic.squad is the better version of those two, while CatmullRom is an intriguing alternative, but may need better treatment of endpoints.

In your equations (A38a) and (A38b) for the quadrangle points, isn't it strange that the contribution by one "chord" is weighted by some time intervals, while the contribution from the other chord is not weighted? To me, this seems very suspicious: changing the chord will change the quadrangle point, but changing the time interval won't?

No, quite the opposite. It might help to point out a notational quirk in my writing. $A_i$ is a control point in the interval between $ti$ and $t{i+1}$, but $Bi$ is a control point in the interval between $t{i-1}$ and $t_i$. (My thinking was that $C(t; R_i, Ai, B{i+1}, t{i+1})$ looks more symmetrical because $B{i+1}$ is hanging out with its friend $t_{i+1}$, and if you're interested in what's happening at $t_i$, you mostly need to worry about $A_i$ and $B_i$.) So the time factor is there to account for the fact that the chord it multiplies is from a different time interval, where the "size" of the chord needs to be rescaled to correct for how long it applies for — since these chords actually represent velocities. When the chord is from the segment in question, we don't need to rescale; when it's from a different segment we do need to rescale.

Actually, I realized that for my random-quaternions example, my method of evaluating the splines on a finer time step was just not up to the task. I've improved that significantly, checked that the result seems to have converged, and updated my gist again. Now, CatmullRom no longer has the lowest torque — that title goes back to quaternionic.squad by a pretty good margin, and Squad is the worst.

To be specific, the average torque of Squad is about 34% larger than that of quaternionic.squad, and the average torque of CatmullRom is about 17% larger.

So I'm fully satisfied with the current state of quaternionic.squad.

Thanks a lot for your detailed response, this is all really interesting!

Sadly, I currently have very little time to look into this, but I want to have a closer look within the next few weeks and then I'll report back.

I still haven't had time to further look into this, but I have very good related news: I have asked at the Computer History Museum if they can scan their copy of the 1987 course notes including Shoemake's paper that introduces Squad. And they promptly did! Since today, the course notes are available as PDF at the bottom of the page https://www.computerhistory.org/collections/catalog/102724883.

Now people can finally read and cite the correct paper!

Ooh, thanks for that!

So, on page 113, he defines the vector $\mathup{t}_n$ to be the average of "the tangents of arcs to adjacent points". I phrased this in terms of the average "velocity" to account for the different time step sizes. To do so, we can divide each logarithm by the time step size relevant to the arc described by that logarithm. So I would rewrite his expression as

$$ \mathup{t}_n = \frac{\mathrm{ln}\left( qn^{-1} q{n+1} \right) / (t{n+1}-t{n}) - \mathrm{ln}\left( qn^{-1} q{n-1} \right) / (t{n}-t{n-1})} {2} $$

Then, (flipping the sign of the second logarithm by taking the inverse inside of it) that turns out to be exactly the fraction on the left-hand sides of my Eqs. (A37). So I would argue that the squad in this package is the most natural extension of Shoemake's original. (Though, again, there is no "right" answer.)

On a related note, if you're still uncomfortable with my Eqs. (A38), it occurs to me that we can rewrite them to look more symmetrical by moving the $\mathbf{R}_i$ to the left-hand side, taking the logarithm of each side, and rearranging the time values so that we're much more clearly relating velocities to each other:

$$ \frac{\log \left( \bar{\mathbf{R}}_i \mathbf{A}_i \right)} {t{i+1}-t{i}} = \frac{1}{4} \left( -\frac{\log \left( \bar{\mathbf{R}}_i \mathbf{R}_{i+1} \right)} {t{i+1}-t{i}} +\frac{\log \left( \bar{\mathbf{R}}_{i-1} \mathbf{R}_{i} \right)} {t{i}-t{i-1}} \right), $$

$$ \frac{\log \left( \bar{\mathbf{R}}_i \mathbf{B}_i \right)} {t{i}-t{i-1}} = \frac{1}{4} \left( -\frac{\log \left( \bar{\mathbf{R}}_{i} \mathbf{R}_{i+1} \right)} {t{i+1}-t{i}} +\frac{\log \left( \bar{\mathbf{R}}_{i-1} \mathbf{R}_{i} \right)} {t{i}-t{i-1}} \right). $$

Note that the right-hand sides are equal, which is what Shoemake shows in the final equation of his Sec. 7; I rewrite his equation as

$$ \frac{\log \left( \bar{\mathbf{R}}_i \mathbf{A}_i \right)} {t{i+1}-t{i}} = \frac{\log \left( \bar{\mathbf{R}}_i \mathbf{B}_i \right)} {t{i}-t{i-1}} = -\frac{ \frac{\log \left( \bar{\mathbf{R}}_{i} \mathbf{R}_{i+1} \right)} {t{i+1}-t{i}} +\frac{\log \left( \bar{\mathbf{R}}_{i} \mathbf{R}_{i-1} \right)} {t{i}-t{i-1}} } {4}. $$

If the time-step sizes are equal, it's easy to see that $\mathbf{A}_i = \mathbf{B}_i$, as he points out.

moble / quaternionic

Equations for Squad? #41