src/main/common/olc.c

   1 #include <math.h>
   2
   3 #include "common/maths.h"
   4 #include "common/olc.h"
   5
   6 // This is a port of https://github.com/google/open-location-code/blob/master/c/olc.c
   7 // to avoid double floating point math and use integer math as much as possible.
   8
   9 #define SEPARATOR_CHAR '+'
  10 #define SEPARATOR_POS 8
  11 #define PADDING_CHAR '0'
  12
  13 #define ENCODING_BASE 20
  14 #define PAIR_CODE_LEN 10u
  15 #define CODE_LEN_MAX 15u
  16
  17 #define GRID_COLS 4
  18 #define GRID_ROWS (ENCODING_BASE / GRID_COLS)
  19
  20
  21 #define LAT_MAX (90 * OLC_DEG_MULTIPLIER)
  22 #define LON_MAX (180 * OLC_DEG_MULTIPLIER)
  23
  24 static const char alphabet[] = "23456789CFGHJMPQRVWX";
  25
  26 // Initialized via init_constants()
  27 static olc_coord_t grid_size;
  28 static olc_coord_t initial_resolution;
  29
  30 static void init_constants(void)
  31 {
  32     static int inited = 0;
  33     if (inited) {
  34         return;
  35     }
  36     inited = 1;
  37
  38     // Work out the encoding base exponent necessary to represent 360 degrees.
  39     olc_coord_t initial_exponent = floorf(logf(2 * (LON_MAX / OLC_DEG_MULTIPLIER)) / logf(ENCODING_BASE));
  40
  41     // Work out the enclosing resolution (in degrees) for the grid algorithm.
  42     grid_size = (1 / powf(ENCODING_BASE, PAIR_CODE_LEN / 2 - (initial_exponent + 1))) * OLC_DEG_MULTIPLIER;
  43
  44     // Work out the initial resolution
  45     initial_resolution = powf(ENCODING_BASE, initial_exponent) * OLC_DEG_MULTIPLIER;
  46 }
  47
  48 // Raises a number to an exponent, handling negative exponents.
  49 static float pow_negf(float base, float exponent)
  50 {
  51     if (exponent == 0) {
  52         return 1;
  53     }
  54
  55     if (exponent > 0) {
  56         return powf(base, exponent);
  57     }
  58
  59     return 1 / powf(base, -exponent);
  60 }
  61
  62 // Compute the latitude precision value for a given code length.  Lengths <= 10
  63 // have the same precision for latitude and longitude, but lengths > 10 have
  64 // different precisions due to the grid method having fewer columns than rows.
  65 static float compute_precision_for_length(int length)
  66 {
  67     // Magic numbers!
  68     if (length <= (int)PAIR_CODE_LEN) {
  69         return pow_negf(ENCODING_BASE, floorf((length / -2) + 2));
  70     }
  71
  72     return pow_negf(ENCODING_BASE, -3) / powf(5, length - (int)PAIR_CODE_LEN);
  73 }
  74
  75 static olc_coord_t adjust_latitude(olc_coord_t lat, size_t code_len)
  76 {
  77     if (lat < -LAT_MAX) {
  78         lat = -LAT_MAX;
  79     }
  80     if (lat > LAT_MAX) {
  81         lat = LAT_MAX;
  82     }
  83     if (lat >= LAT_MAX) {
  84         // Subtract half the code precision to get the latitude into the code area.
  85         olc_coord_t precision = compute_precision_for_length(code_len) * OLC_DEG_MULTIPLIER;
  86         lat -= precision / 2;
  87     }
  88     return lat;
  89 }
  90
  91 static olc_coord_t normalize_longitude(olc_coord_t lon)
  92 {
  93     while (lon < -LON_MAX) {
  94         lon += LON_MAX;
  95         lon += LON_MAX;
  96     }
  97     while (lon >= LON_MAX) {
  98         lon -= LON_MAX;
  99         lon -= LON_MAX;
 100     }
 101     return lon;
 102 }
 103
 104 // Encodes positive range lat,lon into a sequence of OLC lat/lon pairs.  This
 105 // uses pairs of characters (latitude and longitude in that order) to represent
 106 // each step in a 20x20 grid.  Each code, therefore, has 1/400th the area of
 107 // the previous code.
 108 static unsigned encode_pairs(uolc_coord_t lat, uolc_coord_t lon, size_t length, char *buf, size_t bufsize)
 109 {
 110     if ((length + 1) >= bufsize) {
 111         buf[0] = '\0';
 112         return 0;
 113     }
 114
 115     unsigned pos = 0;
 116     olc_coord_t resolution = initial_resolution;
 117     // Add two digits on each pass.
 118     for (size_t digit_count = 0;
 119          digit_count < length;
 120          digit_count += 2, resolution /= ENCODING_BASE) {
 121         size_t digit_value;
 122
 123         // Do the latitude - gets the digit for this place and subtracts that
 124         // for the next digit.
 125         digit_value = lat / resolution;
 126         lat -= digit_value * resolution;
 127         buf[pos++] = alphabet[digit_value];
 128
 129         // Do the longitude - gets the digit for this place and subtracts that
 130         // for the next digit.
 131         digit_value = lon / resolution;
 132         lon -= digit_value * resolution;
 133         buf[pos++] = alphabet[digit_value];
 134
 135         // Should we add a separator here?
 136         if (pos == SEPARATOR_POS && pos < length) {
 137             buf[pos++] = SEPARATOR_CHAR;
 138         }
 139     }
 140     while (pos < SEPARATOR_POS) {
 141         buf[pos++] = PADDING_CHAR;
 142     }
 143     if (pos == SEPARATOR_POS) {
 144         buf[pos++] = SEPARATOR_CHAR;
 145     }
 146     buf[pos] = '\0';
 147     return pos;
 148 }
 149
 150 // Encodes a location using the grid refinement method into an OLC string.  The
 151 // grid refinement method divides the area into a grid of 4x5, and uses a
 152 // single character to refine the area.  The grid squares use the OLC
 153 // characters in order to number the squares as follows:
 154 //
 155 //   R V W X
 156 //   J M P Q
 157 //   C F G H
 158 //   6 7 8 9
 159 //   2 3 4 5
 160 //
 161 // This allows default accuracy OLC codes to be refined with just a single
 162 // character.
 163 static int encode_grid(uolc_coord_t lat, uolc_coord_t lon, size_t length,
 164                        char *buf, size_t bufsize)
 165 {
 166     if ((length + 1) >= bufsize) {
 167         buf[0] = '\0';
 168         return 0;
 169     }
 170
 171     int pos = 0;
 172
 173     olc_coord_t lat_grid_size = grid_size;
 174     olc_coord_t lon_grid_size = grid_size;
 175
 176     lat %= lat_grid_size;
 177     lon %= lon_grid_size;
 178
 179     for (size_t i = 0; i < length; i++) {
 180         olc_coord_t lat_div = lat_grid_size / GRID_ROWS;
 181         olc_coord_t lon_div = lon_grid_size / GRID_COLS;
 182
 183         if (lat_div == 0 || lon_div == 0) {
 184             // This case happens when OLC_DEG_MULTIPLIER doesn't have enough
 185             // precision for the requested length.
 186             break;
 187         }
 188
 189         // Work out the row and column.
 190         size_t row = lat / lat_div;
 191         size_t col = lon / lon_div;
 192         lat_grid_size /= GRID_ROWS;
 193         lon_grid_size /= GRID_COLS;
 194         lat -= row * lat_grid_size;
 195         lon -= col * lon_grid_size;
 196         buf[pos++] = alphabet[row * GRID_COLS + col];
 197     }
 198     buf[pos] = '\0';
 199     return pos;
 200 }
 201
 202 int olc_encode(olc_coord_t lat, olc_coord_t lon, size_t length, char *buf, size_t bufsize)
 203 {
 204     int pos = 0;
 205
 206     length = MIN(length, CODE_LEN_MAX);
 207
 208     // Adjust latitude and longitude so they fall into positive ranges.
 209     uolc_coord_t alat = adjust_latitude(lat, length) + LAT_MAX;
 210     uolc_coord_t alon = normalize_longitude(lon) + LON_MAX;
 211
 212     init_constants();
 213
 214     pos += encode_pairs(alat, alon, MIN(length, PAIR_CODE_LEN), buf + pos, bufsize - pos);
 215     // If the requested length indicates we want grid refined codes.
 216     if (length > PAIR_CODE_LEN) {
 217         pos += encode_grid(alat, alon, length - PAIR_CODE_LEN, buf + pos, bufsize - pos);
 218     }
 219     buf[pos] = '\0';
 220     return pos;
 221 }