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 | */ |
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 | */ |
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 | |
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 | |
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 | |
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 | |
159 | int geohashEncodeWGS84(double longitude, double latitude, uint8_t step, |
160 | GeoHashBits *hash) { |
161 | return geohashEncodeType(longitude, latitude, step, hash); |
162 | } |
163 | |
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 | |
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 | |
202 | int geohashDecodeWGS84(const GeoHashBits hash, GeoHashArea *area) { |
203 | return geohashDecodeType(hash, area); |
204 | } |
205 | |
206 | int 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 | |
217 | int 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 | |
224 | int geohashDecodeToLongLatWGS84(const GeoHashBits hash, double *xy) { |
225 | return geohashDecodeToLongLatType(hash, xy); |
226 | } |
227 | |
228 | static 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 | |
248 | static 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 | |
266 | void 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 | |