How do I apply a digital filter to points on a sphere

Question

Given a set of points on a sphere, how can I implement a higher order low pass filter on them?

At the moment, I am just multiplying the vectors from the input and output set by their weights and summing them, component-wise. While this works well when the angles are small, the filter response isn't right as the angles get larger than 30 degrees. It needs to perfectly compliment the response of a high pass filtered gyroscope, which works with angles, not vectors.

What I'm essentially looking for is an interpolation method for multiple (in my case 8) weighted points on a sphere. The sum of all weights is one. The algorithm needs to be such, that if it were applied to only two vectors $a$ and $b$, with the weights $(1 - t)$ and $(t)$ respectively, a constant rate of change in t would result in a constant angular rate of change in the output. Say, for example we had as our two vectors $a = (1, 0, 0)$ and $b = (0, 1, 0)$, the angle ab is 90 degrees. If we set $t$ to $1/3$, we should expect the result, $c$, to be 30 degrees away from $a$. If we use the standard vector multiplication rules to calculate $c$ as $a*(1 - t) + b*t$, the angle between $a$ and $c$ is closer to 22 degrees. This error gets worse as the angle $ab$ increases. I know that the slerp algorithm exists for just two points, but I need one that works on any number of points.

The problem is that standard linear interpolation between points on a sphere represented as vectors, does not interpolate correctly when you look at the angles.

Answer 1

I once researched the literature on related topics for an internal company report, but I no longer work at that company and I can't remember all the references. One which would at least serve as a good starting point is

Samuel R. Buss and Jay Fillmore
Spherical Averages and Applications to Spherical Splines and Interpolation
ACM Transactions on Graphics 20 (2001) 95-126.

In fact, skim-reading it I recognise that their spherical average is the one we ended up using, so you can consider this a recommendation.

Answer 2

I don't think you've actually specified the behaviour you want for more than two vectors. In that case, there's no obvious way to represent the vectors in terms of angles between them. One way to attempt to generalize your specification for two vectors would be to require that the interpolation works as if it were done two vectors at a time. Unfortunately, that doesn't work because the result would depend on the order in which you use the vectors.

For instance, take the three vectors $(1,0,0)$, $(0,1,0)$, $(0,0,1)$ and combine them with weight $1/3$ each. By symmetry, the result must be $(1,1,1)/\sqrt3$. If we do this two vectors at a time, we can first combine $(1,0,0)$ and $(0,1,0)$ with equal weight to get $(1,1,0)/\sqrt2$. For this to work, we should then get $(1,1,1)/\sqrt3$ by combining $(1,1,0)/\sqrt2$ and $(0,0,1)$ with weights $2/3$ and $1/3$, respectively. We don't however, as you can see by calculating the cosines of $(1,1,1)/\sqrt3$ with those vectors, which are $\sqrt{2/3}$ and $\sqrt{1/3}$, respectively, corresponding to angles of about $35.3^\circ$ and $54.7^\circ$, respectively, and not $30^\circ$ and $60^\circ$ as they would have to be.

So this isn't a well-defined problem so far, and I don't see how to naturally turn it into one.

[Edit:]

I'm wondering from your comment whether perhaps you're only interested in vectors that all lie on a great circle and don't care what the algorithm does otherwise. In that case, you could find the vector most nearly orthogonal to all the vectors, which is the eigenvector of the coordinate covariance matrix $\sum_ix_{i\alpha}x_{i\beta}$ with least eigenvalue, rotate that into the $z$ axis, find the polar coordinates $\theta_i,\phi_i$ of the vectors in the rotated system, take the weighted average of the $\theta$ values, do whatever you want with the $\phi$ values as if they were angles on a circle, and rotate the result back. If your vectors all lie on a great circle, they'll get rotated onto the equator and their $\phi$ angles will behave as you want them to.

Note that this will fail in the above example, since those three vectors would have $\phi$ values $2\pi/3$ apart, and there would be no sensible way to average those; so this is really only a solution for the case where the vectors approximately lie on a great circle and are concentrated in one half of it.