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 (65536).
50 * From: https://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN
51 */
interleave64(uint32_t xlo,uint32_t ylo)52 static 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 */
deinterleave64(uint64_t interleaved)82 static 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
geohashGetCoordRange(GeoHashRange * long_range,GeoHashRange * lat_range)112 void 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
geohashEncode(const GeoHashRange * long_range,const GeoHashRange * lat_range,double longitude,double latitude,uint8_t step,GeoHashBits * hash)121 int 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
geohashEncodeType(double longitude,double latitude,uint8_t step,GeoHashBits * hash)153 int 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
geohashEncodeWGS84(double longitude,double latitude,uint8_t step,GeoHashBits * hash)159 int geohashEncodeWGS84(double longitude, double latitude, uint8_t step,
160 GeoHashBits *hash) {
161 return geohashEncodeType(longitude, latitude, step, hash);
162 }
163
geohashDecode(const GeoHashRange long_range,const GeoHashRange lat_range,const GeoHashBits hash,GeoHashArea * area)164 int 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
geohashDecodeType(const GeoHashBits hash,GeoHashArea * area)196 int 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
geohashDecodeWGS84(const GeoHashBits hash,GeoHashArea * area)202 int geohashDecodeWGS84(const GeoHashBits hash, GeoHashArea *area) {
203 return geohashDecodeType(hash, area);
204 }
205
geohashDecodeAreaToLongLat(const GeoHashArea * area,double * xy)206 int geohashDecodeAreaToLongLat(const GeoHashArea *area, double *xy) {
207 if (!xy) return 0;
208 xy[0] = (area->longitude.min + area->longitude.max) / 2;
209 xy[1] = (area->latitude.min + area->latitude.max) / 2;
210 return 1;
211 }
212
geohashDecodeToLongLatType(const GeoHashBits hash,double * xy)213 int geohashDecodeToLongLatType(const GeoHashBits hash, double *xy) {
214 GeoHashArea area = {{0}};
215 if (!xy || !geohashDecodeType(hash, &area))
216 return 0;
217 return geohashDecodeAreaToLongLat(&area, xy);
218 }
219
geohashDecodeToLongLatWGS84(const GeoHashBits hash,double * xy)220 int geohashDecodeToLongLatWGS84(const GeoHashBits hash, double *xy) {
221 return geohashDecodeToLongLatType(hash, xy);
222 }
223
geohash_move_x(GeoHashBits * hash,int8_t d)224 static void geohash_move_x(GeoHashBits *hash, int8_t d) {
225 if (d == 0)
226 return;
227
228 uint64_t x = hash->bits & 0xaaaaaaaaaaaaaaaaULL;
229 uint64_t y = hash->bits & 0x5555555555555555ULL;
230
231 uint64_t zz = 0x5555555555555555ULL >> (64 - hash->step * 2);
232
233 if (d > 0) {
234 x = x + (zz + 1);
235 } else {
236 x = x | zz;
237 x = x - (zz + 1);
238 }
239
240 x &= (0xaaaaaaaaaaaaaaaaULL >> (64 - hash->step * 2));
241 hash->bits = (x | y);
242 }
243
geohash_move_y(GeoHashBits * hash,int8_t d)244 static void geohash_move_y(GeoHashBits *hash, int8_t d) {
245 if (d == 0)
246 return;
247
248 uint64_t x = hash->bits & 0xaaaaaaaaaaaaaaaaULL;
249 uint64_t y = hash->bits & 0x5555555555555555ULL;
250
251 uint64_t zz = 0xaaaaaaaaaaaaaaaaULL >> (64 - hash->step * 2);
252 if (d > 0) {
253 y = y + (zz + 1);
254 } else {
255 y = y | zz;
256 y = y - (zz + 1);
257 }
258 y &= (0x5555555555555555ULL >> (64 - hash->step * 2));
259 hash->bits = (x | y);
260 }
261
geohashNeighbors(const GeoHashBits * hash,GeoHashNeighbors * neighbors)262 void geohashNeighbors(const GeoHashBits *hash, GeoHashNeighbors *neighbors) {
263 neighbors->east = *hash;
264 neighbors->west = *hash;
265 neighbors->north = *hash;
266 neighbors->south = *hash;
267 neighbors->south_east = *hash;
268 neighbors->south_west = *hash;
269 neighbors->north_east = *hash;
270 neighbors->north_west = *hash;
271
272 geohash_move_x(&neighbors->east, 1);
273 geohash_move_y(&neighbors->east, 0);
274
275 geohash_move_x(&neighbors->west, -1);
276 geohash_move_y(&neighbors->west, 0);
277
278 geohash_move_x(&neighbors->south, 0);
279 geohash_move_y(&neighbors->south, -1);
280
281 geohash_move_x(&neighbors->north, 0);
282 geohash_move_y(&neighbors->north, 1);
283
284 geohash_move_x(&neighbors->north_west, -1);
285 geohash_move_y(&neighbors->north_west, 1);
286
287 geohash_move_x(&neighbors->north_east, 1);
288 geohash_move_y(&neighbors->north_east, 1);
289
290 geohash_move_x(&neighbors->south_east, 1);
291 geohash_move_y(&neighbors->south_east, -1);
292
293 geohash_move_x(&neighbors->south_west, -1);
294 geohash_move_y(&neighbors->south_west, -1);
295 }
296