Init commit.HEADmaster

1 files changed, 219 insertions, 0 deletions

diff --git a/src/geohash_helper.c b/src/geohash_helper.c
new file mode 100644
index 0000000..77d8ab3
--- /dev/null
+++ b/src/geohash_helper.c
@@ -0,0 +1,219 @@
+/*
+ * Copyright (c) 2013-2014, yinqiwen <[email protected]>
+ * Copyright (c) 2014, Matt Stancliff <[email protected]>.
+ * Copyright (c) 2015-2016, Salvatore Sanfilippo <[email protected]>.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ *  * Redistributions of source code must retain the above copyright notice,
+ *    this list of conditions and the following disclaimer.
+ *  * Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *  * Neither the name of Redis nor the names of its contributors may be used
+ *    to endorse or promote products derived from this software without
+ *    specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+ * THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/* This is a C++ to C conversion from the ardb project.
+ * This file started out as:
+ * https://github.com/yinqiwen/ardb/blob/d42503/src/geo/geohash_helper.cpp
+ */
+
+#include "fmacros.h"
+#include "geohash_helper.h"
+#include "debugmacro.h"
+#include <math.h>
+
+#define D_R (M_PI / 180.0)
+#define R_MAJOR 6378137.0
+#define R_MINOR 6356752.3142
+#define RATIO (R_MINOR / R_MAJOR)
+#define ECCENT (sqrt(1.0 - (RATIO *RATIO)))
+#define COM (0.5 * ECCENT)
+
+/// @brief The usual PI/180 constant
+const double DEG_TO_RAD = 0.017453292519943295769236907684886;
+/// @brief Earth's quatratic mean radius for WGS-84
+const double EARTH_RADIUS_IN_METERS = 6372797.560856;
+
+const double MERCATOR_MAX = 20037726.37;
+const double MERCATOR_MIN = -20037726.37;
+
+static inline double deg_rad(double ang) { return ang * D_R; }
+static inline double rad_deg(double ang) { return ang / D_R; }
+
+/* This function is used in order to estimate the step (bits precision)
+ * of the 9 search area boxes during radius queries. */
+uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
+    if (range_meters == 0) return 26;
+    int step = 1;
+    while (range_meters < MERCATOR_MAX) {
+        range_meters *= 2;
+        step++;
+    }
+    step -= 2; /* Make sure range is included in most of the base cases. */
+
+    /* Wider range torwards the poles... Note: it is possible to do better
+     * than this approximation by computing the distance between meridians
+     * at this latitude, but this does the trick for now. */
+    if (lat > 66 || lat < -66) {
+        step--;
+        if (lat > 80 || lat < -80) step--;
+    }
+
+    /* Frame to valid range. */
+    if (step < 1) step = 1;
+    if (step > 26) step = 26;
+    return step;
+}
+
+/* Return the bounding box of the search area centered at latitude,longitude
+ * having a radius of radius_meter. bounds[0] - bounds[2] is the minimum
+ * and maxium longitude, while bounds[1] - bounds[3] is the minimum and
+ * maximum latitude. */
+int geohashBoundingBox(double longitude, double latitude, double radius_meters,
+                       double *bounds) {
+    if (!bounds) return 0;
+
+    bounds[0] = longitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
+    bounds[2] = longitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
+    bounds[1] = latitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
+    bounds[3] = latitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
+    return 1;
+}
+
+/* Return a set of areas (center + 8) that are able to cover a range query
+ * for the specified position and radius. */
+GeoHashRadius geohashGetAreasByRadius(double longitude, double latitude, double radius_meters) {
+    GeoHashRange long_range, lat_range;
+    GeoHashRadius radius;
+    GeoHashBits hash;
+    GeoHashNeighbors neighbors;
+    GeoHashArea area;
+    double min_lon, max_lon, min_lat, max_lat;
+    double bounds[4];
+    int steps;
+
+    geohashBoundingBox(longitude, latitude, radius_meters, bounds);
+    min_lon = bounds[0];
+    min_lat = bounds[1];
+    max_lon = bounds[2];
+    max_lat = bounds[3];
+
+    steps = geohashEstimateStepsByRadius(radius_meters,latitude);
+
+    geohashGetCoordRange(&long_range,&lat_range);
+    geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
+    geohashNeighbors(&hash,&neighbors);
+    geohashDecode(long_range,lat_range,hash,&area);
+
+    /* Check if the step is enough at the limits of the covered area.
+     * Sometimes when the search area is near an edge of the
+     * area, the estimated step is not small enough, since one of the
+     * north / south / west / east square is too near to the search area
+     * to cover everything. */
+    int decrease_step = 0;
+    {
+        GeoHashArea north, south, east, west;
+
+        geohashDecode(long_range, lat_range, neighbors.north, &north);
+        geohashDecode(long_range, lat_range, neighbors.south, &south);
+        geohashDecode(long_range, lat_range, neighbors.east, &east);
+        geohashDecode(long_range, lat_range, neighbors.west, &west);
+
+        if (geohashGetDistance(longitude,latitude,longitude,north.latitude.max)
+            < radius_meters) decrease_step = 1;
+        if (geohashGetDistance(longitude,latitude,longitude,south.latitude.min)
+            < radius_meters) decrease_step = 1;
+        if (geohashGetDistance(longitude,latitude,east.longitude.max,latitude)
+            < radius_meters) decrease_step = 1;
+        if (geohashGetDistance(longitude,latitude,west.longitude.min,latitude)
+            < radius_meters) decrease_step = 1;
+    }
+
+    if (steps > 1 && decrease_step) {
+        steps--;
+        geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
+        geohashNeighbors(&hash,&neighbors);
+        geohashDecode(long_range,lat_range,hash,&area);
+    }
+
+    /* Exclude the search areas that are useless. */
+    if (area.latitude.min < min_lat) {
+        GZERO(neighbors.south);
+        GZERO(neighbors.south_west);
+        GZERO(neighbors.south_east);
+    }
+    if (area.latitude.max > max_lat) {
+        GZERO(neighbors.north);
+        GZERO(neighbors.north_east);
+        GZERO(neighbors.north_west);
+    }
+    if (area.longitude.min < min_lon) {
+        GZERO(neighbors.west);
+        GZERO(neighbors.south_west);
+        GZERO(neighbors.north_west);
+    }
+    if (area.longitude.max > max_lon) {
+        GZERO(neighbors.east);
+        GZERO(neighbors.south_east);
+        GZERO(neighbors.north_east);
+    }
+    radius.hash = hash;
+    radius.neighbors = neighbors;
+    radius.area = area;
+    return radius;
+}
+
+GeoHashRadius geohashGetAreasByRadiusWGS84(double longitude, double latitude,
+                                           double radius_meters) {
+    return geohashGetAreasByRadius(longitude, latitude, radius_meters);
+}
+
+GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash) {
+    uint64_t bits = hash.bits;
+    bits <<= (52 - hash.step * 2);
+    return bits;
+}
+
+/* Calculate distance using haversin great circle distance formula. */
+double geohashGetDistance(double lon1d, double lat1d, double lon2d, double lat2d) {
+    double lat1r, lon1r, lat2r, lon2r, u, v;
+    lat1r = deg_rad(lat1d);
+    lon1r = deg_rad(lon1d);
+    lat2r = deg_rad(lat2d);
+    lon2r = deg_rad(lon2d);
+    u = sin((lat2r - lat1r) / 2);
+    v = sin((lon2r - lon1r) / 2);
+    return 2.0 * EARTH_RADIUS_IN_METERS *
+           asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
+}
+
+int geohashGetDistanceIfInRadius(double x1, double y1,
+                                 double x2, double y2, double radius,
+                                 double *distance) {
+    *distance = geohashGetDistance(x1, y1, x2, y2);
+    if (*distance > radius) return 0;
+    return 1;
+}
+
+int geohashGetDistanceIfInRadiusWGS84(double x1, double y1, double x2,
+                                      double y2, double radius,
+                                      double *distance) {
+    return geohashGetDistanceIfInRadius(x1, y1, x2, y2, radius, distance);
+}