1/*
2 * Copyright (c) 2013-2014, yinqiwen <[email protected]>
3 * Copyright (c) 2014, Matt Stancliff <[email protected]>.
4 * Copyright (c) 2015-2016, Salvatore Sanfilippo <[email protected]>.
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions are met:
9 *
10 * * Redistributions of source code must retain the above copyright notice,
11 * this list of conditions and the following disclaimer.
12 * * Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * * Neither the name of Redis nor the names of its contributors may be used
16 * to endorse or promote products derived from this software without
17 * specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
23 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
29 * THE POSSIBILITY OF SUCH DAMAGE.
30 */
31#include "geohash.h"
32
33/**
34 * Hashing works like this:
35 * Divide the world into 4 buckets. Label each one as such:
36 * -----------------
37 * | | |
38 * | | |
39 * | 0,1 | 1,1 |
40 * -----------------
41 * | | |
42 * | | |
43 * | 0,0 | 1,0 |
44 * -----------------
45 */
46
47/* Interleave lower bits of x and y, so the bits of x
48 * are in the even positions and bits from y in the odd;
49 * x and y must initially be less than 2**32 (4294967296).
50 * From: https://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN
51 */
52static inline uint64_t interleave64(uint32_t xlo, uint32_t ylo) {
53 static const uint64_t B[] = {0x5555555555555555ULL, 0x3333333333333333ULL,
54 0x0F0F0F0F0F0F0F0FULL, 0x00FF00FF00FF00FFULL,
55 0x0000FFFF0000FFFFULL};
56 static const unsigned int S[] = {1, 2, 4, 8, 16};
57
58 uint64_t x = xlo;
59 uint64_t y = ylo;
60
61 x = (x | (x << S[4])) & B[4];
62 y = (y | (y << S[4])) & B[4];
63
64 x = (x | (x << S[3])) & B[3];
65 y = (y | (y << S[3])) & B[3];
66
67 x = (x | (x << S[2])) & B[2];
68 y = (y | (y << S[2])) & B[2];
69
70 x = (x | (x << S[1])) & B[1];
71 y = (y | (y << S[1])) & B[1];
72
73 x = (x | (x << S[0])) & B[0];
74 y = (y | (y << S[0])) & B[0];
75
76 return x | (y << 1);
77}
78
79/* reverse the interleave process
80 * derived from http://stackoverflow.com/questions/4909263
81 */
82static inline uint64_t deinterleave64(uint64_t interleaved) {
83 static const uint64_t B[] = {0x5555555555555555ULL, 0x3333333333333333ULL,
84 0x0F0F0F0F0F0F0F0FULL, 0x00FF00FF00FF00FFULL,
85 0x0000FFFF0000FFFFULL, 0x00000000FFFFFFFFULL};
86 static const unsigned int S[] = {0, 1, 2, 4, 8, 16};
87
88 uint64_t x = interleaved;
89 uint64_t y = interleaved >> 1;
90
91 x = (x | (x >> S[0])) & B[0];
92 y = (y | (y >> S[0])) & B[0];
93
94 x = (x | (x >> S[1])) & B[1];
95 y = (y | (y >> S[1])) & B[1];
96
97 x = (x | (x >> S[2])) & B[2];
98 y = (y | (y >> S[2])) & B[2];
99
100 x = (x | (x >> S[3])) & B[3];
101 y = (y | (y >> S[3])) & B[3];
102
103 x = (x | (x >> S[4])) & B[4];
104 y = (y | (y >> S[4])) & B[4];
105
106 x = (x | (x >> S[5])) & B[5];
107 y = (y | (y >> S[5])) & B[5];
108
109 return x | (y << 32);
110}
111
112void geohashGetCoordRange(GeoHashRange *long_range, GeoHashRange *lat_range) {
113 /* These are constraints from EPSG:900913 / EPSG:3785 / OSGEO:41001 */
114 /* We can't geocode at the north/south pole. */
115 long_range->max = GEO_LONG_MAX;
116 long_range->min = GEO_LONG_MIN;
117 lat_range->max = GEO_LAT_MAX;
118 lat_range->min = GEO_LAT_MIN;
119}
120
121int geohashEncode(const GeoHashRange *long_range, const GeoHashRange *lat_range,
122 double longitude, double latitude, uint8_t step,
123 GeoHashBits *hash) {
124 /* Check basic arguments sanity. */
125 if (hash == NULL || step > 32 || step == 0 ||
126 RANGEPISZERO(lat_range) || RANGEPISZERO(long_range)) return 0;
127
128 /* Return an error when trying to index outside the supported
129 * constraints. */
130 if (longitude > GEO_LONG_MAX || longitude < GEO_LONG_MIN ||
131 latitude > GEO_LAT_MAX || latitude < GEO_LAT_MIN) return 0;
132
133 hash->bits = 0;
134 hash->step = step;
135
136 if (latitude < lat_range->min || latitude > lat_range->max ||
137 longitude < long_range->min || longitude > long_range->max) {
138 return 0;
139 }
140
141 double lat_offset =
142 (latitude - lat_range->min) / (lat_range->max - lat_range->min);
143 double long_offset =
144 (longitude - long_range->min) / (long_range->max - long_range->min);
145
146 /* convert to fixed point based on the step size */
147 lat_offset *= (1ULL << step);
148 long_offset *= (1ULL << step);
149 hash->bits = interleave64(lat_offset, long_offset);
150 return 1;
151}
152
153int geohashEncodeType(double longitude, double latitude, uint8_t step, GeoHashBits *hash) {
154 GeoHashRange r[2] = {{0}};
155 geohashGetCoordRange(&r[0], &r[1]);
156 return geohashEncode(&r[0], &r[1], longitude, latitude, step, hash);
157}
158
159int geohashEncodeWGS84(double longitude, double latitude, uint8_t step,
160 GeoHashBits *hash) {
161 return geohashEncodeType(longitude, latitude, step, hash);
162}
163
164int geohashDecode(const GeoHashRange long_range, const GeoHashRange lat_range,
165 const GeoHashBits hash, GeoHashArea *area) {
166 if (HASHISZERO(hash) || NULL == area || RANGEISZERO(lat_range) ||
167 RANGEISZERO(long_range)) {
168 return 0;
169 }
170
171 area->hash = hash;
172 uint8_t step = hash.step;
173 uint64_t hash_sep = deinterleave64(hash.bits); /* hash = [LAT][LONG] */
174
175 double lat_scale = lat_range.max - lat_range.min;
176 double long_scale = long_range.max - long_range.min;
177
178 uint32_t ilato = hash_sep; /* get lat part of deinterleaved hash */
179 uint32_t ilono = hash_sep >> 32; /* shift over to get long part of hash */
180
181 /* divide by 2**step.
182 * Then, for 0-1 coordinate, multiply times scale and add
183 to the min to get the absolute coordinate. */
184 area->latitude.min =
185 lat_range.min + (ilato * 1.0 / (1ull << step)) * lat_scale;
186 area->latitude.max =
187 lat_range.min + ((ilato + 1) * 1.0 / (1ull << step)) * lat_scale;
188 area->longitude.min =
189 long_range.min + (ilono * 1.0 / (1ull << step)) * long_scale;
190 area->longitude.max =
191 long_range.min + ((ilono + 1) * 1.0 / (1ull << step)) * long_scale;
192
193 return 1;
194}
195
196int geohashDecodeType(const GeoHashBits hash, GeoHashArea *area) {
197 GeoHashRange r[2] = {{0}};
198 geohashGetCoordRange(&r[0], &r[1]);
199 return geohashDecode(r[0], r[1], hash, area);
200}
201
202int geohashDecodeWGS84(const GeoHashBits hash, GeoHashArea *area) {
203 return geohashDecodeType(hash, area);
204}
205
206int geohashDecodeAreaToLongLat(const GeoHashArea *area, double *xy) {
207 if (!xy) return 0;
208 xy[0] = (area->longitude.min + area->longitude.max) / 2;
209 if (xy[0] > GEO_LONG_MAX) xy[0] = GEO_LONG_MAX;
210 if (xy[0] < GEO_LONG_MIN) xy[0] = GEO_LONG_MIN;
211 xy[1] = (area->latitude.min + area->latitude.max) / 2;
212 if (xy[1] > GEO_LAT_MAX) xy[1] = GEO_LAT_MAX;
213 if (xy[1] < GEO_LAT_MIN) xy[1] = GEO_LAT_MIN;
214 return 1;
215}
216
217int geohashDecodeToLongLatType(const GeoHashBits hash, double *xy) {
218 GeoHashArea area = {{0}};
219 if (!xy || !geohashDecodeType(hash, &area))
220 return 0;
221 return geohashDecodeAreaToLongLat(&area, xy);
222}
223
224int geohashDecodeToLongLatWGS84(const GeoHashBits hash, double *xy) {
225 return geohashDecodeToLongLatType(hash, xy);
226}
227
228static void geohash_move_x(GeoHashBits *hash, int8_t d) {
229 if (d == 0)
230 return;
231
232 uint64_t x = hash->bits & 0xaaaaaaaaaaaaaaaaULL;
233 uint64_t y = hash->bits & 0x5555555555555555ULL;
234
235 uint64_t zz = 0x5555555555555555ULL >> (64 - hash->step * 2);
236
237 if (d > 0) {
238 x = x + (zz + 1);
239 } else {
240 x = x | zz;
241 x = x - (zz + 1);
242 }
243
244 x &= (0xaaaaaaaaaaaaaaaaULL >> (64 - hash->step * 2));
245 hash->bits = (x | y);
246}
247
248static void geohash_move_y(GeoHashBits *hash, int8_t d) {
249 if (d == 0)
250 return;
251
252 uint64_t x = hash->bits & 0xaaaaaaaaaaaaaaaaULL;
253 uint64_t y = hash->bits & 0x5555555555555555ULL;
254
255 uint64_t zz = 0xaaaaaaaaaaaaaaaaULL >> (64 - hash->step * 2);
256 if (d > 0) {
257 y = y + (zz + 1);
258 } else {
259 y = y | zz;
260 y = y - (zz + 1);
261 }
262 y &= (0x5555555555555555ULL >> (64 - hash->step * 2));
263 hash->bits = (x | y);
264}
265
266void geohashNeighbors(const GeoHashBits *hash, GeoHashNeighbors *neighbors) {
267 neighbors->east = *hash;
268 neighbors->west = *hash;
269 neighbors->north = *hash;
270 neighbors->south = *hash;
271 neighbors->south_east = *hash;
272 neighbors->south_west = *hash;
273 neighbors->north_east = *hash;
274 neighbors->north_west = *hash;
275
276 geohash_move_x(&neighbors->east, 1);
277 geohash_move_y(&neighbors->east, 0);
278
279 geohash_move_x(&neighbors->west, -1);
280 geohash_move_y(&neighbors->west, 0);
281
282 geohash_move_x(&neighbors->south, 0);
283 geohash_move_y(&neighbors->south, -1);
284
285 geohash_move_x(&neighbors->north, 0);
286 geohash_move_y(&neighbors->north, 1);
287
288 geohash_move_x(&neighbors->north_west, -1);
289 geohash_move_y(&neighbors->north_west, 1);
290
291 geohash_move_x(&neighbors->north_east, 1);
292 geohash_move_y(&neighbors->north_east, 1);
293
294 geohash_move_x(&neighbors->south_east, 1);
295 geohash_move_y(&neighbors->south_east, -1);
296
297 geohash_move_x(&neighbors->south_west, -1);
298 geohash_move_y(&neighbors->south_west, -1);
299}
300