Help save net neutrality! Learn more.
Close

Weighted Spherical Mean Calculations

Help
Anonymous
2013-06-16
2013-10-19
  • Anonymous - 2013-06-16

    Hello all,

    Great library so thank you!
    Have a question, trying to calculate the weighted spherical mean of points and was wondering how I could achieve this with this library. To be exact, I am trying to do something similar to this: http://research.microsoft.com/pubs/136532/volley.pdf on page 7.

    Thank you in advance!
    Nick

     
  • Charles Karney

    Charles Karney - 2013-06-17

    I recommend the following book/paper on how to do weighted spherical means:

    • K. V. Mardia & P. E. Jupp, "Directional Statistics", (Wiley, 1999).
    • S. R. Buss & J. P. Fillmore, "Spherical Averages and Applications to Spherical Splines and Interpolation", ACM Transactions on Graphics 20, 95-126 (2001).

    The method used in the Volley paper is what Buss & Fillmore call a "progressive slerp". There are a couple of problems with this method:

    1. The implementation given in the Volley paper is horrible. The first equation for d is well known to be ill-conditioned for nearby points. Then I got an atan(0/0) in the equation for gamma with phi_A = 90, phi_B = 0 (note that the authors use phi for the colatitude). Finally, with w = 0.5, I don't get the same answer if I interchange x_A and x_B. The problems are easily fixed using GeographicLib (and you get the ellipsoidal generalization with no extra trouble); see the function slerpmean below (and I've replaced the recursion with a simple for loop).
    2. The bigger problem with progressive slerp is that the result you get depends on the order of the points. This is nearly always undesirable.

    There are two other methods which I can recommend which don't depend on the order of the points.

    1. Euclidean mean; this is discussed in Mardia & Jupp. It's easy to implement and easy to generalize to an ellipsoid. It allows big problems to be split into pieces and the results combined. The only criticism I know of is that the result depends on how the ellipsoid surface is embedded in 3d. This is unlikely to be a concern for geodesic applications. This is implemented by euclidmean below. A by product of this method is an estimate of the variance of the points (via h).
    2. Geodesic mean as given by Buss & Fillmore. This minimizes the RMS geodesic distance from the points to the mean. This is an intrinsic method (does not depend on the embedding). However, it's iterative and it's not possible to solve by combining the results of separate pieces. The function geodmean below is an ellipsoidal generalization of Buss & Fillmore's Algorithm A1.

    The three methods are illustrated by this code. It reads in a sequence of points (lat lon) (assumed here to be all of equal weight) and prints out the three means. I recommend you use euclidmean.

    #include <iostream>
    #include <iomanip>
    #include <vector>
    #include <GeographicLib/Geodesic.hpp>
    #include <GeographicLib/Geocentric.hpp>
    #include <GeographicLib/AzimuthalEquidistant.hpp>
    
    using namespace GeographicLib;
    using namespace std;
    
    // Hold postion + weight
    struct pt {
      double lat, lon, w;
      pt(double latx, double lonx, double wx = 1)
        : lat(latx)
        , lon(lonx)
        , w(wx) {}
    };
    
    // Weighted mean of two points via slerp (generalized to the ellipsoid)
    pt slerp(const Geodesic& g, const pt& p1, const pt& p2) {
      double s12, azi1, azi2;
      g.Inverse(p1.lat, p1.lon, p2.lat, p2.lon, s12, azi1, azi2);
      double lat, lon;
      g.Direct(p1.lat, p1.lon, azi1, s12 * p2.w/(p1.w + p2.w), lat, lon);
      return pt(lat, lon, p1.w + p2.w);
    }
    
    // Mean of several points via progressive slerp
    pt slerpmean(const Geodesic& g, const vector<pt> points) {
      pt p(0, 0, 0);
      for (size_t i = 0; i < points.size(); ++i)
        p = slerp(g, points[i], p);
      return p;
    }
    
    // Euclidean mean of points.
    // returned value of h (height of mean point) is measure of variance.
    pt euclidmean(const Geocentric& c, const vector<pt> points, double& h) {
      double X = 0, Y = 0, Z = 0, W = 0;
      for (size_t i = 0; i < points.size(); ++i) {
        double x, y, z, w;
        c.Forward(points[i].lat, points[i].lon, 0, x, y, z);
        w = points[i].w;
        X += w * x; Y += w * y; Z += w * z; W += w;
      }
      X /= W; Y /= W; Z /= W;
      double lat, lon;
      c.Reverse(X, Y, Z, lat, lon, h);
      return pt(lat, lon, W);
    }
    
    // Mean of points miminizing RMS geodesic distance to points.
    // lat0,lon0 is starting guess.
    // See Buss and Fillmore, "Spherical Averages and Applications to Spherical
    // Splines and Interpolation", ACM Transactions on Graphics 20 (2001) 95-126.
    pt geodmean(const Geodesic& g, const vector<pt> points,
                double lat0, double lon0) {
      const AzimuthalEquidistant a(g);
      double W = 0;
      for (size_t k = 0; k < 20; ++k) {
        double X = 0, Y = 0;
        W = 0;
        for (size_t i = 0; i < points.size(); ++i) {
          double x, y, w;
          a.Forward(lat0, lon0, points[i].lat, points[i].lon, x, y);
          w = points[i].w;
          X += w * x; Y += w * y; W += w;
        }
        X /= W; Y /= W;
        a.Reverse(lat0, lon0, X, Y, lat0, lon0);
      }
      return pt(lat0, lon0, W);
    }
    
    int main() {
      double a = Constants::WGS84_a(), f = Constants::WGS84_f();
      const Geodesic g(a, f);
      const Geocentric c(a, f);
      vector<pt> points;
      double lat, lon;
      while (cin >> lat >> lon)
        points.push_back(pt(lat,lon));
      cout << fixed;
    
      pt sp = slerpmean(g, points);
      cout << setprecision(8) << sp.lat << " " << sp.lon << "\n";
    
      double h;
      pt ep = euclidmean(c, points, h);
      cout << setprecision(8) << ep.lat << " " << ep.lon << " "
           << setprecision(3) << h << "\n";
    
      pt gp = geodmean(g, points, ep.lat, ep.lon);
      cout << setprecision(8) << gp.lat << " " << gp.lon << "\n";
    }
    
     
  • Anonymous - 2013-06-17

    Wow thank you for the great explanation!!

    That most definitely helps with the implementation of a Volley like system - only trying to estimate it for a comparison. I was aware that their explanation and overview of what they do was a little off and with a poor choice for an algorithm, so good to see someone else agree.

    I did some digging and also know of these two algorithms. Will go with your recommendation, the Mardia & Jupp algo.

    Again thank you very much for your time and help on this!

     
  • Anonymous - 2013-07-25

    Quick question - looking into the citations you gave I couldn't find the complexity of these two algorithms. Also trying to port this code you supplied into the Java version of your library - any pointers?

    Many thanks!

     
  • Charles Karney

    Charles Karney - 2013-07-26

    Euclidean mean is O(N) where N is the number of points. Geodesic mean
    is O(N) times the number of iterations. I don't have a tight bound on
    this. However for reasonable restricted sets of points (all within
    10000 km, say) the convergence appears to be quadratic (like Newton's
    method) so the number of iterations is on the order of the logarithm of
    the number of bits of precision (typically 5 iterations are needed).

    I should also mention that the geodesic mean is older than the Buss &
    Fillmore paper suggests. Look at the Wikipedia article on "Frechet
    mean" (or "Karcher mean").

     
  • Anonymous - 2013-08-14

    Thank you Charles, that was very helpful. Wanted to make sure this wasn't going to kill my performance when dealing with relatively large numbers. Just adding this into a system I've built now and noticed Geocentric has not been implemented for the Java library. is that correct? if so then I believe I can only use the worst of the three; slerpmean.

     
  • Charles Karney

    Charles Karney - 2013-08-14

    I recommend you implement GeoCentric in java. The code for the pretty short and you might simplify you job by skipping the option of returning the rotation matrices. Then you are left with implementing two pretty short functions IntForward and IntReverse.

     


Anonymous

Cancel  Add attachments