xapian-core/geospatial/latlong_metrics.cc

   1 /** @file
   2  * @brief Geospatial distance metrics.
   3  */
   4 /* Copyright 2008 Lemur Consulting Ltd
   5  * Copyright 2011 Richard Boulton
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License as
   9  * published by the Free Software Foundation; either version 2 of the
  10  * License, or (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
  20  * USA
  21  */
  22
  23 #include <config.h>
  24
  25 #include "xapian/geospatial.h"
  26 #include "xapian/error.h"
  27 #include "serialise-double.h"
  28
  29 #include <cmath>
  30
  31 using namespace Xapian;
  32 using namespace std;
  33
  34 /** Quadratic mean radius of the Earth in metres.
  35  */
  36 #define QUAD_EARTH_RADIUS_METRES 6372797.6
  37
  38 LatLongMetric::~LatLongMetric()
  39 {
  40 }
  41
  42 double
  43 LatLongMetric::operator()(const LatLongCoords& a,
  44                           const LatLongCoords& b) const
  45 {
  46     if (a.empty() || b.empty()) {
  47         throw InvalidArgumentError("Empty coordinate list supplied to LatLongMetric::operator()()");
  48     }
  49     double min_dist = 0.0;
  50     bool have_min = false;
  51     for (LatLongCoordsIterator a_iter = a.begin();
  52          a_iter != a.end();
  53          ++a_iter)
  54     {
  55         for (LatLongCoordsIterator b_iter = b.begin();
  56              b_iter != b.end();
  57              ++b_iter)
  58         {
  59             double dist = pointwise_distance(*a_iter, *b_iter);
  60             if (!have_min) {
  61                 min_dist = dist;
  62                 have_min = true;
  63             } else if (dist < min_dist) {
  64                 min_dist = dist;
  65             }
  66         }
  67     }
  68     return min_dist;
  69 }
  70
  71 double
  72 LatLongMetric::operator()(const LatLongCoords& a,
  73                           const char* b_ptr, size_t b_len) const
  74 {
  75     if (a.empty() || b_len == 0) {
  76         throw InvalidArgumentError("Empty coordinate list supplied to LatLongMetric::operator()()");
  77     }
  78     double min_dist = 0.0;
  79     bool have_min = false;
  80     LatLongCoord b;
  81     const char * b_end = b_ptr + b_len;
  82     while (b_ptr != b_end) {
  83         b.unserialise(&b_ptr, b_end);
  84         for (LatLongCoordsIterator a_iter = a.begin();
  85              a_iter != a.end();
  86              ++a_iter)
  87         {
  88             double dist = pointwise_distance(*a_iter, b);
  89             if (!have_min) {
  90                 min_dist = dist;
  91                 have_min = true;
  92             } else if (dist < min_dist) {
  93                 min_dist = dist;
  94             }
  95         }
  96     }
  97     return min_dist;
  98 }
  99
 100
 101 GreatCircleMetric::GreatCircleMetric()
 102         : radius(QUAD_EARTH_RADIUS_METRES)
 103 {}
 104
 105 GreatCircleMetric::GreatCircleMetric(double radius_)
 106         : radius(radius_)
 107 {}
 108
 109 double
 110 GreatCircleMetric::pointwise_distance(const LatLongCoord& a,
 111                                       const LatLongCoord& b) const
 112 {
 113     double lata = a.latitude * (M_PI / 180.0);
 114     double latb = b.latitude * (M_PI / 180.0);
 115
 116     double latdiff = lata - latb;
 117     double longdiff = (a.longitude - b.longitude) * (M_PI / 180.0);
 118
 119     double sin_half_lat = sin(latdiff / 2);
 120     double sin_half_long = sin(longdiff / 2);
 121     double h = sin_half_lat * sin_half_lat +
 122             sin_half_long * sin_half_long * cos(lata) * cos(latb);
 123     if (rare(h > 1.0)) {
 124         // Clamp to 1.0, asin(1.0) = M_PI / 2.0.
 125         return radius * M_PI;
 126     }
 127     return 2 * radius * asin(sqrt(h));
 128 }
 129
 130 LatLongMetric *
 131 GreatCircleMetric::clone() const
 132 {
 133     return new GreatCircleMetric(radius);
 134 }
 135
 136 string
 137 GreatCircleMetric::name() const
 138 {
 139     return "Xapian::GreatCircleMetric";
 140 }
 141
 142 string
 143 GreatCircleMetric::serialise() const
 144 {
 145     return serialise_double(radius);
 146 }
 147
 148 LatLongMetric *
 149 GreatCircleMetric::unserialise(const string& s) const
 150 {
 151     const char * p = s.data();
 152     const char * end = p + s.size();
 153
 154     double new_radius = unserialise_double(&p, end);
 155     if (p != end) {
 156         throw Xapian::NetworkError("Bad serialised GreatCircleMetric - junk at end");
 157     }
 158
 159     return new GreatCircleMetric(new_radius);
 160 }