Skip to content

Commit

Permalink
Merge faster encode/decode implementation from @Steve132
Browse files Browse the repository at this point in the history
  • Loading branch information
@yinqiwen
yinqiwen committed Sep 16, 2014
1 parent bfddf70 commit 8840c4c
Show file tree
Hide file tree
Showing 2 changed files with 205 additions and 71 deletions.
262 changes: 195 additions & 67 deletions geohash.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
*/
#include <stddef.h>
#include <stdio.h>
#include <math.h>
#include "geohash.h"

/**
Expand All @@ -44,16 +45,13 @@
* -----------------
*/

int geohash_encode(const GeoHashRange* lat_r, const GeoHashRange* lon_r, double latitude, double longitude,
uint8_t step, GeoHashBits* hash)
int geohash_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitude, double longitude, uint8_t step,
GeoHashBits* hash)
{
if (NULL == hash || step > 32 || step == 0 || NULL == lat_r || NULL == lon_r)
if (NULL == hash || step > 32 || step == 0)
{
return -1;
}
GeoHashRange lat_range, lon_range;
lat_range = *lat_r;
lon_range = *lon_r;
hash->bits = 0;
hash->step = step;
uint8_t i = 0;
Expand Down Expand Up @@ -98,24 +96,23 @@ static inline uint8_t get_bit(uint64_t bits, uint8_t pos)
return (bits >> pos) & 0x01;
}

int geohash_decode(const GeoHashRange* lat_range, const GeoHashRange* lon_range, const GeoHashBits* hash,
GeoHashArea* area)
int geohash_decode(GeoHashRange lat_range, GeoHashRange lon_range, GeoHashBits hash, GeoHashArea* area)
{
if (NULL == hash || NULL == area || NULL == lat_range || NULL == lon_range)
if (NULL == area)
{
return -1;
}
area->hash = *hash;
area->hash = hash;
uint8_t i = 0;
area->latitude.min = lat_range->min;
area->latitude.max = lat_range->max;
area->longitude.min = lon_range->min;
area->longitude.max = lon_range->max;
for (; i < hash->step; i++)
area->latitude.min = lat_range.min;
area->latitude.max = lat_range.max;
area->longitude.min = lon_range.min;
area->longitude.max = lon_range.max;
for (; i < hash.step; i++)
{
uint8_t lat_bit, lon_bit;
lon_bit = get_bit(hash->bits, (hash->step - i) * 2 - 1);
lat_bit = get_bit(hash->bits, (hash->step - i) * 2 - 2);
lon_bit = get_bit(hash.bits, (hash.step - i) * 2 - 1);
lat_bit = get_bit(hash.bits, (hash.step - i) * 2 - 2);
if (lat_bit == 0)
{
area->latitude.max = (area->latitude.max + area->latitude.min) / 2;
Expand All @@ -136,6 +133,132 @@ int geohash_decode(const GeoHashRange* lat_range, const GeoHashRange* lon_range,
return 0;
}

static inline uint64_t interleave64(uint32_t xlo, uint32_t ylo)
{
static const uint64_t B[] =
{ 0x5555555555555555, 0x3333333333333333, 0x0F0F0F0F0F0F0F0F, 0x00FF00FF00FF00FF, 0x0000FFFF0000FFFF };
static const unsigned int S[] =
{ 1, 2, 4, 8, 16 };

uint64_t x = xlo; // Interleave lower bits of x and y, so the bits of x
uint64_t y = ylo; // are in the even positions and bits from y in the odd; //https://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN

// x and y must initially be less than 2**32.

x = (x | (x << S[4])) & B[4];
y = (y | (y << S[4])) & B[4];

x = (x | (x << S[3])) & B[3];
y = (y | (y << S[3])) & B[3];

x = (x | (x << S[2])) & B[2];
y = (y | (y << S[2])) & B[2];

x = (x | (x << S[1])) & B[1];
y = (y | (y << S[1])) & B[1];

x = (x | (x << S[0])) & B[0];
y = (y | (y << S[0])) & B[0];

return x | (y << 1);
}

static inline uint64_t deinterleave64(uint64_t interleaved)
{
static const uint64_t B[] =
{ 0x5555555555555555, 0x3333333333333333, 0x0F0F0F0F0F0F0F0F, 0x00FF00FF00FF00FF, 0x0000FFFF0000FFFF,
0x00000000FFFFFFFF };
static const unsigned int S[] =
{ 0, 1, 2, 4, 8, 16 };

uint64_t x = interleaved; ///reverse the interleave process (http://stackoverflow.com/questions/4909263/how-to-efficiently-de-interleave-bits-inverse-morton)
uint64_t y = interleaved >> 1;

x = (x | (x >> S[0])) & B[0];
y = (y | (y >> S[0])) & B[0];

x = (x | (x >> S[1])) & B[1];
y = (y | (y >> S[1])) & B[1];

x = (x | (x >> S[2])) & B[2];
y = (y | (y >> S[2])) & B[2];

x = (x | (x >> S[3])) & B[3];
y = (y | (y >> S[3])) & B[3];

x = (x | (x >> S[4])) & B[4];
y = (y | (y >> S[4])) & B[4];

x = (x | (x >> S[5])) & B[5];
y = (y | (y >> S[5])) & B[5];

return x | (y << 32);
}


int geohash_fast_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitude, double longitude, uint8_t step,
GeoHashBits* hash)
{
if (NULL == hash || step > 32 || step == 0)
{
return -1;
}
hash->bits = 0;
hash->step = step;
uint8_t i = 0;
if (latitude < lat_range.min || latitude > lat_range.max || longitude < lon_range.min || longitude > lon_range.max)
{
return -1;
}

//the algorithm computes the morton code for the geohash location within the range.
//this can be done MUCH more efficiently using the following code

//compute the coordinate in the range 0-1
double lat_offset = (latitude - lat_range.min) / (lat_range.max - lat_range.min);
double lon_offset = (longitude - lon_range.min) / (lon_range.max - lon_range.min);

//convert it to fixed point based on the step size
lat_offset *= (1 << step);
lon_offset *= (1 << step);

uint32_t ilato = (uint32_t) lat_offset;
uint32_t ilono = (uint32_t) lon_offset;

//interleave the bits to create the morton code. No branching and no bounding
hash->bits = interleave64(ilato, ilono);
return 0;
}
int geohash_fast_decode(GeoHashRange lat_range, GeoHashRange lon_range, GeoHashBits hash, GeoHashArea* area)
{
if (NULL == area)
{
return -1;
}
area->hash = hash;
uint8_t step = hash.step;
uint64_t xyhilo = deinterleave64(hash.bits); //decode morton code

double lat_scale = lat_range.max - lat_range.min;
double lon_scale = lon_range.max - lon_range.min;

uint32_t ilato = xyhilo; //get back the original integer coordinates
uint32_t ilono = xyhilo >> 32;

//double lat_offset=ilato;
//double lon_offset=ilono;
//lat_offset /= (1<<step);
//lon_offset /= (1<<step);

//the ldexp call converts the integer to a double,then divides by 2**step to get the 0-1 coordinate, which is then multiplied times scale and added to the min to get the absolute coordinate
area->latitude.min = lat_range.min + ldexp(ilato, -step) * lat_scale;
area->latitude.max = lat_range.min + ldexp(ilato + 1, -step) * lat_scale;
area->longitude.min = lon_range.min + ldexp(ilono, -step) * lon_scale;
area->longitude.max = lon_range.min + ldexp(ilono + 1, -step) * lon_scale;

return 0;
}

static int geohash_move_x(GeoHashBits* hash, int8_t d)
{
if (d == 0)
Expand Down Expand Up @@ -180,7 +303,7 @@ static int geohash_move_y(GeoHashBits* hash, int8_t d)
return 0;
}

int geohash_get_neighbors(const GeoHashBits* hash, GeoHashNeighbors* neighbors)
int geohash_get_neighbors(GeoHashBits hash, GeoHashNeighbors* neighbors)
{
geohash_get_neighbor(hash, GEOHASH_NORTH, &neighbors->north);
geohash_get_neighbor(hash, GEOHASH_EAST, &neighbors->east);
Expand All @@ -193,13 +316,13 @@ int geohash_get_neighbors(const GeoHashBits* hash, GeoHashNeighbors* neighbors)
return 0;
}

int geohash_get_neighbor(const GeoHashBits* hash, GeoDirection direction, GeoHashBits* neighbor)
int geohash_get_neighbor(GeoHashBits hash, GeoDirection direction, GeoHashBits* neighbor)
{
if (NULL == neighbor)
{
return -1;
}
*neighbor = *hash;
*neighbor = hash;
switch (direction)
{
case GEOHASH_NORTH:
Expand Down Expand Up @@ -292,51 +415,56 @@ GeoHashBits geohash_next_righttop(GeoHashBits bits)
}

/*
int main()
{
GeoHashBits hash;
GeoHashNeighbors neighbors;
GeoHashRange lat_range, lon_range;
lat_range.max = 20037726.37;
lat_range.min = -20037726.37;
lon_range.max = 20037726.37;
lon_range.min = -20037726.37;
double radius = 5000;
double latitude = 9741705.20;
double longitude = 5417390.90;
GeoHashBits hash_min, hash_max, hash_lt, hash_gr;
geohash_encode(&lat_range, &lon_range, latitude, longitude, 13, &hash);
geohash_encode(&lat_range, &lon_range, latitude - radius, longitude - radius, 13, &hash_min);
geohash_encode(&lat_range, &lon_range, latitude + radius, longitude + radius, 13, &hash_max);
geohash_encode(&lat_range, &lon_range, latitude + radius, longitude - radius, 13, &hash_lt);
geohash_encode(&lat_range, &lon_range, latitude - radius, longitude + radius, 13, &hash_gr);
printf("%lld\n", hash.bits);
geohash_get_neighbors(&hash, &neighbors);
printf("%lld\n", hash.bits);
printf("%lld\n", neighbors.east.bits);
printf("%lld\n", neighbors.west.bits);
printf("%lld\n", neighbors.south.bits);
printf("%lld\n", neighbors.north.bits);
printf("%lld\n", neighbors.north_west.bits);
printf("%lld\n", neighbors.north_east.bits);
printf("%lld\n", neighbors.south_east.bits);
printf("%lld\n", neighbors.south_west.bits);
printf("##%lld\n", hash_min.bits);
printf("##%lld\n", hash_max.bits);
printf("##%lld\n", hash_lt.bits);
printf("##%lld\n", hash_gr.bits);
// geohash_encode(&lat_range, &lon_range, 9741705.20, 5417390.90, 13, &hash);
// printf("from %lld to %lld \n", hash.bits << 2, (hash.bits+1) << 2);
GeoHashArea area;
geohash_decode(&lat_range, &lon_range, &hash, &area);
printf("%.10f %.10f\n", area.latitude.min, area.latitude.max);
printf("%.10f %.10f\n", area.longitude.min, area.longitude.max);
return 0;
}
int main()
{
GeoHashBits hash, fast_hash;
GeoHashNeighbors neighbors;
GeoHashRange lat_range, lon_range;
lat_range.max = 20037726.37;
lat_range.min = -20037726.37;
lon_range.max = 20037726.37;
lon_range.min = -20037726.37;
double radius = 5000;
double latitude = 9741705.20;
double longitude = 5417390.90;
GeoHashBits hash_min, hash_max, hash_lt, hash_gr;
geohash_encode(lat_range, lon_range, latitude, longitude, 24, &hash);
geohash_fast_encode(lat_range, lon_range, latitude, longitude, 24, &fast_hash);
geohash_encode(lat_range, lon_range, latitude - radius, longitude - radius, 13, &hash_min);
geohash_encode(lat_range, lon_range, latitude + radius, longitude + radius, 13, &hash_max);
geohash_encode(lat_range, lon_range, latitude + radius, longitude - radius, 13, &hash_lt);
geohash_encode(lat_range, lon_range, latitude - radius, longitude + radius, 13, &hash_gr);
printf("## %lld\n", hash.bits);
printf("## %lld\n", fast_hash.bits);
geohash_get_neighbors(hash, &neighbors);
printf("%lld\n", hash.bits);
printf("%lld\n", neighbors.east.bits);
printf("%lld\n", neighbors.west.bits);
printf("%lld\n", neighbors.south.bits);
printf("%lld\n", neighbors.north.bits);
printf("%lld\n", neighbors.north_west.bits);
printf("%lld\n", neighbors.north_east.bits);
printf("%lld\n", neighbors.south_east.bits);
printf("%lld\n", neighbors.south_west.bits);
printf("##%lld\n", hash_min.bits);
printf("##%lld\n", hash_max.bits);
printf("##%lld\n", hash_lt.bits);
printf("##%lld\n", hash_gr.bits);
// geohash_encode(&lat_range, &lon_range, 9741705.20, 5417390.90, 13, &hash);
// printf("from %lld to %lld \n", hash.bits << 2, (hash.bits+1) << 2);
GeoHashArea area, area1;
geohash_decode(lat_range, lon_range, hash, &area);
geohash_fast_decode(lat_range, lon_range, hash, &area1);
printf("%.10f %.10f\n", area.latitude.min, area.latitude.max);
printf("%.10f %.10f\n", area.longitude.min, area.longitude.max);
printf("%.10f %.10f\n", area1.latitude.min, area1.latitude.max);
printf("%.10f %.10f\n", area1.longitude.min, area1.longitude.max);
return 0;
}
*/
14 changes: 10 additions & 4 deletions geohash.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,11 +84,17 @@ extern "C"
* 0:success
* -1:failed
*/
int geohash_encode(const GeoHashRange* lat_range, const GeoHashRange* lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash);
int geohash_decode(const GeoHashRange* lat_range, const GeoHashRange* lon_range, const GeoHashBits* hash, GeoHashArea* area);
int geohash_get_neighbors(const GeoHashBits* hash, GeoHashNeighbors* neighbors);
int geohash_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash);
int geohash_decode(GeoHashRange lat_range, GeoHashRange lon_range, GeoHashBits hash, GeoHashArea* area);

int geohash_get_neighbor(const GeoHashBits* hash, GeoDirection direction, GeoHashBits* neighbor);
/*
* Fast encode/decode version, more magic in implementation.
*/
int geohash_fast_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash);
int geohash_fast_decode(GeoHashRange lat_range, GeoHashRange lon_range, GeoHashBits hash, GeoHashArea* area);

int geohash_get_neighbors(GeoHashBits hash, GeoHashNeighbors* neighbors);
int geohash_get_neighbor(GeoHashBits hash, GeoDirection direction, GeoHashBits* neighbor);

GeoHashBits geohash_next_leftbottom(GeoHashBits bits);
GeoHashBits geohash_next_rightbottom(GeoHashBits bits);
Expand Down

0 comments on commit 8840c4c

Please sign in to comment.