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

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.

-
Can you explain how your points are assigned velocities and weights, in your final remark? Is this purely a problem of points on a sphere, or are there other aspects? – Joseph O'Rourke Oct 2 '11 at 22:48
I've clarified a bit, I hope it helps. The problem is simply that I cannot find a method of combining weighted points in spherical geometry. When the angles are small, it can be approximated using Euclidean geometry like I'm doing at the moment. – Hannesh Oct 3 '11 at 9:17

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.

-
@Hannesh: This looks good; I'd use this instead, as it handles all cases gracefully. I wouldn't mind if you switch the accepted answer to this :-) – joriki Oct 3 '11 at 12:33
This is exactly what I needed! And thanks for being so modest joriki ;) – Hannesh Oct 3 '11 at 13:17

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.

-
I see your point, I tried your example myself with the slerp algorithm too and didn't get the intended result either. What I want as the behaviour is that if we applied the filter to a set of vectors such as { (1, 0, 0), (1, 0, 0), (0, 1, 0), (0, 1, 0) ... }, the response of the angles should be exactly the same as if we filtered the set { 0°, 0°, 90°, 90° ... }. I know this is defined on the surface of a circle in 2D, but are you insisting that it is just not possible in spherical geometry? – Hannesh Oct 3 '11 at 11:25
No, this is certainly possible; what I'm saying is that this doesn't specify how the filter should operate on general vectors. For instance, in my example, this specifies what the result should be if one of the weights is zero (i.e. on the edges of the spherical triangle formed by the vectors), but it says nothing about what the result should be if all three weights are non-zero. There are many possible continuations of the result on the edges to the interior of the triangle, and I don't see a particularly natural one. Or perhaps you don't care as long it works for vectors on a great circle? – joriki Oct 3 '11 at 11:35
@Hannesh: I've added a solution to the answer that will work if you don't care what happens when the vectors aren't on a great circle. – joriki Oct 3 '11 at 11:57
True, I'm having trouble imagining what should happen in other cases. That could be good enough, I'll give that one a go. Thanks! - Also, it's not that I don't care what the algorithm does otherwise, it should still give reasonable results. I'll see how it goes – Hannesh Oct 3 '11 at 12:02