Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Slightly improved API. Dramatically improved implementation #1

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 103 additions & 67 deletions geohash.c
Original file line number Diff line number Diff line change
Expand Up @@ -43,53 +43,96 @@
* | 0,0 | 1,0 |
* -----------------
*/
//updated files.
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;

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_encode(const GeoHashRange* lat_r, const GeoHashRange* lon_r, double latitude, double longitude,
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;
if (latitude < lat_range.min || latitude > lat_range.max || longitude < lon_range.min || longitude > lon_range.max)
{
return -1;
}

for (; i < step; i++)
{
uint8_t lat_bit, lon_bit;
if (lat_range.max - latitude >= latitude - lat_range.min)
{
lat_bit = 0;
lat_range.max = (lat_range.max + lat_range.min) / 2;
}
else
{
lat_bit = 1;
lat_range.min = (lat_range.max + lat_range.min) / 2;
}
if (lon_range.max - longitude >= longitude - lon_range.min)
{
lon_bit = 0;
lon_range.max = (lon_range.max + lon_range.min) / 2;
}
else
{
lon_bit = 1;
lon_range.min = (lon_range.max + lon_range.min) / 2;
}
hash->bits <<= 1;
hash->bits += lon_bit;
hash->bits <<= 1;
hash->bits += lat_bit;
}


//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;
}

Expand All @@ -98,41 +141,34 @@ 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,
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;
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++)
{
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);
if (lat_bit == 0)
{
area->latitude.max = (area->latitude.max + area->latitude.min) / 2;
}
else
{
area->latitude.min = (area->latitude.max + area->latitude.min) / 2;
}
if (lon_bit == 0)
{
area->longitude.max = (area->longitude.max + area->longitude.min) / 2;
}
else
{
area->longitude.min = (area->longitude.max + area->longitude.min) / 2;
}
}
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;

ilato=xyhilo; //get back the original integer coordinates
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;
}

Expand Down Expand Up @@ -180,7 +216,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 +229,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
8 changes: 4 additions & 4 deletions geohash.h
Original file line number Diff line number Diff line change
Expand Up @@ -84,11 +84,11 @@ 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_neighbors(GeoHashBits hash, GeoHashNeighbors* neighbors);

int geohash_get_neighbor(const GeoHashBits* hash, GeoDirection direction, GeoHashBits* neighbor);
int geohash_get_neighbor(GeoHashBits hash, GeoDirection direction, GeoHashBits* neighbor);

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