geometry.c

/* Copyright 2013 Bliksem Labs. See the LICENSE file at the top-level directory of this distribution and at https://github.com/bliksemlabs/rrrr/. */

/* geometry.c */
#include "geometry.h"

#include "config.h"
#include <math.h>
#include <stdio.h>

// Mean of Earth's equatorial and meridional circumferences.
#define EARTH_CIRCUMFERENCE 40041438.5
// UINT32_MAX is also the full range of INT32.
#define INT32_RANGE UINT32_MAX
// We could have more resolution in the latitude direction by mapping 90 degrees to the int32 range
// instead of 180, but keeping both axes at the same scale enables efficent distance calculations.
// In any case the extra Y resolution is unnecessary, since 1 brad is already just under 1cm.
#define METERS_PER_BRAD (EARTH_CIRCUMFERENCE / INT32_RANGE)
// Must be scaled according to latitude for use in the longitude direction.
#define METERS_PER_DEGREE_LAT (EARTH_CIRCUMFERENCE / 360.0)

static inline double radians (double degrees) {
    return degrees * M_PI / 180;
}

static inline double degrees (double radians) {
    return radians * 180 / M_PI;
}

// Sinusoidal / equirectangular projection.
// Replace with fast estimation polynomial.
static inline double xscale_at_y (uint32_t y_brads) {
    return cos(y_brads * M_PI * 2 / UINT32_MAX);
}

static inline double xscale_at_lat (double latitude) {
    return cos(radians(latitude));
}

static inline double coord_diff_meters (int32_t o1, int32_t o2) {
     return (o2 - o1) * METERS_PER_BRAD;
}

/*
Several mappings from lat/lon to integers to are possible.
The three major decisions are:
1. Signed or unsigned brads.
2. Zero uint32_t in brads is at the prime meridian or its antipode.
3. Whether or not to store scaled x values.

Design decisions and reflections:
If we map zero brads to 0 longitude and xvalues are not scaled, unsigned and signed brads have identical binary representations.
The main advantage of using unsigned brads is that is prevents hashgrids from "reflecting" around the origin when normal division is performed.
The unsigned conversion could be performed within the hashgrid itself, and binning could be performed with masks (shift out bits within cell, zero out high bits, rest is cell number).
If the x values are stored pre-scaled we will not get proper wrapping at the date line, but bins will be of a consistent size.
If the x values are stored unscaled then we will get seamless wrapping at the date line, but bins will not be of a consistent size.
Not all architectures will use twos complement signed integers... which don't?
In order to efficiently calculate distances, we need the x and y coordinates to be in the same units / at the same scale.
We can either sacrifice vertical resolution and store both coordinates at the x scale, or perform an additional multiply operation at every distance calculation to restore the y coordinate to the same scale.
2. Either we start x=0 at lon=-180 or lon=0.
*/
void coord_from_lat_lon (coord_t *coord, double lat, double lon) {
    coord->y = lat * UINT32_MAX / 360.0;
    coord->x = lon * UINT32_MAX / 360.0 * xscale_at_lat (lat);
}

void coord_from_meters (coord_t *coord, double meters_x, double meters_y) {
    coord->x = meters_x / METERS_PER_BRAD;
    coord->y = meters_y / METERS_PER_BRAD;
}

void coord_from_latlon (coord_t *coord, latlon_t *latlon) {
    coord_from_lat_lon (coord, latlon->lat, latlon->lon);
}

// Returns a monotone increasing function of the distance between points that is faster to calculate.
// This is an order-preserving map, so it can be used to establish closer/farther relationships.
// It is implemented as the square of the distance in latitude-scaled/projected binary radians.
double coord_ersatz_distance (coord_t *c1, coord_t *c2) {
    double dx = c2->x - c1->x;
    double dy = c2->y - c1->y;
    return dx * dx + dy * dy;
} // rename to coord_distance_ersatz, add ersatz_from_meters and meters_from_ersatz

// Returns the equivalent ersatz distance for a given number of meters, for distance comparisons.
double ersatz_distance (double meters) {
    double d_brads = meters / METERS_PER_BRAD ;
    return d_brads * d_brads;
}

double coord_distance_meters (coord_t *c1, coord_t *c2) {
    double dxm = coord_diff_meters(c1->x, c2->x);
    double dym = coord_diff_meters(c1->y, c2->y);
    return sqrt(dxm * dxm + dym * dym);
}

double latlon_distance_meters (latlon_t *ll1, latlon_t *ll2) {
    // Rather than finding and using the average latitude for longitude scaling, or scaling
    // the two latlons separately, just convert the to the internal projected representation.
    coord_t c1, c2;
    coord_from_latlon (&c1, ll1);
    coord_from_latlon (&c2, ll2);
    return coord_distance_meters (&c1, &c2);
}

void latlon_from_coord (latlon_t *latlon, coord_t *coord) {
    latlon->lat = coord->y * 180.0 / INT32_MAX ;
    latlon->lon = coord->x * 180.0 / INT32_MAX / xscale_at_y (coord->y);
}

void latlon_dump (latlon_t *latlon) {
    printf("latlon lat=%f lon=%f \n", latlon->lat, latlon->lon);
}

void coord_dump (coord_t *coord) {
    printf("coordinate x=%d y=%d \n", coord->x, coord->y);
}