Skip to content

Commit c68a13b

Browse files
mapcodefoundationmapcodefoundation
committed
2.1.0 Fix floating point error accumulation
Prevent encode(decode(M)) != M
1 parent 606bc87 commit c68a13b

File tree

2 files changed

+196
-72
lines changed

2 files changed

+196
-72
lines changed

‎mapcodelib/mapcoder.c‎

Lines changed: 192 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -64,8 +64,13 @@ typedef struct {
6464
int context; // input territory context (or negative)
6565
const char *iso; // input territory alphacode (context)
6666
// output
67-
point result; // result
68-
point32 coord32; // result in integer arithmetic (millionts of degrees)
67+
point result; // result
68+
point32 coord32; // result in integer arithmetic (millionts of degrees)
69+
#ifdef FORCE_RECODE
70+
#define MICROMETER (0.000001) // millionth millionth degree ~ 0,1 micron ~ 1/810000th
71+
point range_min;
72+
point range_max;
73+
#endif
6974
} decodeRec;
7075

7176

@@ -301,8 +306,14 @@ static void encodeExtension(char *result, int extrax4, int extray, int dividerx4
301306
const encodeRec *enc) // append extra characters to result for more precision
302307
{
303308
char *s = result + strlen(result);
304-
double encx = (extrax4 + 4 * enc->fraclon) / (dividerx4);
305-
double ency = (extray + enc->fraclat * ydirection) / (dividery);
309+
int factorx = MAX_PRECISION_FACTOR * dividerx4; // 30^4
310+
int factory = MAX_PRECISION_FACTOR * dividery; // 30^4
311+
int valx = (int) floor(MAX_PRECISION_FACTOR * (extrax4 + 4 * enc->fraclon));
312+
int valy = (int) floor(MAX_PRECISION_FACTOR * (extray + enc->fraclat * ydirection));
313+
314+
// protect against floating point errors
315+
if (valx<0) { valx=0; } else if (valx>=factorx) { valx=factorx-1; }
316+
if (valy<0) { valy=0; } else if (valy>=factory) { valy=factory-1; }
306317

307318
if (extraDigits < 0) { extraDigits = 0; } else if (extraDigits > MAX_PRECISION_DIGITS) {
308319
extraDigits = MAX_PRECISION_DIGITS;
@@ -315,20 +326,14 @@ static void encodeExtension(char *result, int extrax4, int extray, int dividerx4
315326
while (extraDigits-- > 0) {
316327
int gx, gy, column1, column2, row1, row2;
317328

318-
encx *= 30;
319-
gx = (int) encx;
320-
if (gx < 0) {
321-
gx = 0;
322-
} else if (gx > 29) {
323-
gx = 29;
324-
}
325-
ency *= 30;
326-
gy = (int) ency;
327-
if (gy < 0) {
328-
gy = 0;
329-
} else if (gy > 29) {
330-
gy = 29;
331-
}
329+
factorx /= 30;
330+
gx = (valx / factorx);
331+
valx -= factorx * gx;
332+
333+
factory /= 30;
334+
gy = (valy / factory);
335+
valy -= factory * gy;
336+
332337
column1 = (gx / 6);
333338
column2 = (gx % 6);
334339
row1 = (gy / 5);
@@ -339,22 +344,21 @@ static void encodeExtension(char *result, int extrax4, int extray, int dividerx4
339344
*s++ = encode_chars[row2 * 6 + column2];
340345
}
341346
*s = 0;
342-
encx -= gx;
343-
ency -= gy;
344347
}
345348
}
346349

347350
#define decodeChar(c) decode_chars[(unsigned char)c] // force c to be in range of the index, between 0 and 255
348351

349352
// this routine takes the integer-arithmeteic decoding results (in millionths of degrees), adds any floating-point precision digits, and returns the result (still in millionths)
350-
static int decodeExtension(decodeRec *dec, int dividerx4, int dividery, int ydirection) {
353+
static int decodeExtension(decodeRec *dec, int dividerx4, int dividery0, int lon_offset4) {
351354
const char *extrapostfix = dec->extension;
352-
double dividerx = dividerx4 / 4.0, processor = 1.0;
353-
dec->result.lon = 0;
354-
dec->result.lat = 0;
355+
double dividerx = dividerx4 / 4.0, dividery = dividery0;
356+
int lon32 = 0;
357+
int lat32 = 0;
358+
int processor = 1;
359+
int odd = 0;
355360
while (*extrapostfix) {
356361
int column1, row1, column2, row2;
357-
double halfcolumn = 0;
358362
int c1 = *extrapostfix++;
359363
c1 = decodeChar(c1);
360364
if (c1 < 0 || c1 == 30) { return -1; } // illegal extension character
@@ -367,25 +371,39 @@ static int decodeExtension(decodeRec *dec, int dividerx4, int dividery, int ydir
367371
column2 = (c2 % 6);
368372
}
369373
else {
370-
row2 = 2;
371-
halfcolumn = 0.5;
372-
column2 = 3;
374+
odd = 1;
375+
row2 = 0;
376+
column2 = 0;
373377
}
374378

375-
processor *= 30;
376-
dec->result.lon += ((column1 * 6 + column2)) / processor;
377-
dec->result.lat += ((row1 * 5 + row2 - halfcolumn)) / processor;
379+
processor *= 30;
380+
lon32 = lon32 * 30 + column1 * 6 + column2;
381+
lat32 = lat32 * 30 + row1 * 5 + row2;
378382
}
379383

380-
dec->result.lon += (0.5 / processor);
381-
dec->result.lat += (0.5 / processor);
382-
383-
dec->result.lon *= dividerx;
384-
dec->result.lat *= (dividery * ydirection);
384+
dec->result.lon = dec->coord32.lon + (lon32 * dividerx / processor) + lon_offset4 / 4.0;
385+
dec->result.lat = dec->coord32.lat + (lat32 * dividery / processor);
385386

386-
dec->result.lon += dec->coord32.lon;
387-
dec->result.lat += dec->coord32.lat;
387+
#ifdef FORCE_RECODE
388+
dec->range_min.lon = dec->result.lon;
389+
dec->range_min.lat = dec->result.lat;
390+
if (odd) {
391+
dec->range_max.lon = dec->range_min.lon + (dividerx / (processor / 6));
392+
dec->range_max.lat = dec->range_min.lat + (dividery / (processor / 5));
393+
} else {
394+
dec->range_max.lon = dec->range_min.lon + (dividerx / processor);
395+
dec->range_max.lat = dec->range_min.lat + (dividery / processor);
396+
} //
397+
#endif // FORCE_RECODE
388398

399+
if (odd) {
400+
dec->result.lon += dividerx / (2 * processor / 6);
401+
dec->result.lat += dividery / (2 * processor / 5);
402+
} else {
403+
dec->result.lon += (dividerx / (2 * processor));
404+
dec->result.lat += (dividery / (2 * processor));
405+
} // not odd
406+
389407
// also convert back to int
390408
dec->coord32.lon = (int) floor(dec->result.lon);
391409
dec->coord32.lat = (int) floor(dec->result.lat);
@@ -597,12 +615,29 @@ static int decodeGrid(decodeRec *dec, int m, int hasHeaderLetter) {
597615

598616
dec->coord32.lon = relx + (difx * dividerx);
599617
dec->coord32.lat = rely + (dify * dividery);
618+
if (!fitsInside(&dec->coord32,m)) {
619+
return -912;
620+
}
600621

601622
{
602-
int err = decodeExtension(dec, dividerx << 2, dividery, 1); // grid
623+
int err = decodeExtension(dec, dividerx << 2, dividery, 0); // grid
603624
if (err) {
604625
return err;
605626
}
627+
#ifdef FORCE_RECODE
628+
if (dec->result.lon >= relx + xgridsize) {
629+
dec->coord32.lon = (int) floor(dec->result.lon = relx + xgridsize - MICROMETER);
630+
} // keep in inner cell
631+
if (dec->result.lat >= rely + ygridsize) {
632+
dec->coord32.lat = (int) floor(dec->result.lat = rely + ygridsize - MICROMETER);
633+
} // keep in inner cell
634+
if (dec->result.lon >= b->maxx) {
635+
dec->coord32.lon = (int) floor(dec->result.lon = b->maxx - MICROMETER);
636+
} // keep in territory
637+
if (dec->result.lat >= b->maxy) {
638+
dec->coord32.lat = (int) floor(dec->result.lat = b->maxy - MICROMETER);
639+
} // keep in territory
640+
#endif // FORCE_RECODE
606641

607642
return 0;
608643
}
@@ -875,7 +910,7 @@ static int decodeNameless(decodeRec *dec, int m) {
875910

876911

877912
if (dx >= xSIDE) {
878-
return -123; //EVER?
913+
return -123;
879914
}
880915

881916
{
@@ -887,8 +922,17 @@ static int decodeNameless(decodeRec *dec, int m) {
887922
dec->coord32.lon = b->minx + ((dx * dividerx4) / 4);
888923
dec->coord32.lat = b->maxy - (dy * dividery);
889924

890-
err = decodeExtension(dec, dividerx4, dividery, -1); // nameless
891-
dec->result.lon += ((dx * dividerx4) % 4) / 4.0;
925+
err = decodeExtension(dec, dividerx4, -dividery, ((dx * dividerx4) % 4)); // nameless
926+
927+
#ifdef FORCE_RECODE
928+
// keep within outer rect
929+
if (dec->result.lat < b->miny) {
930+
dec->result.lat = (dec->coord32.lat = b->miny);
931+
} // keep in territory
932+
if (dec->result.lon >= b->maxx) {
933+
dec->coord32.lon = (int) floor(dec->result.lon = b->maxx - MICROMETER);
934+
} // keep in territory
935+
#endif
892936

893937
return err;
894938
}
@@ -1157,11 +1201,20 @@ static int decodeAutoHeader(decodeRec *dec, int m) {
11571201
dec->coord32.lat > b->maxy) // *** CAREFUL! do this test BEFORE adding remainder...
11581202
{
11591203
return -122; // invalid code
1160-
}
1204+
}
11611205
}
11621206
}
11631207

1164-
err = decodeExtension(dec, dividerx << 2, dividery, -1); // autoheader decode
1208+
err = decodeExtension(dec, dividerx << 2, -dividery, 0); // autoheader decode
1209+
1210+
#ifdef FORCE_RECODE
1211+
if (dec->result.lat < b->miny) {
1212+
dec->result.lat = (dec->coord32.lat = b->miny);
1213+
} // keep in territory
1214+
if (dec->result.lon >= b->maxx) {
1215+
dec->coord32.lon = (int) floor(dec->result.lon = b->maxx - MICROMETER);
1216+
} // keep in territory
1217+
#endif
11651218

11661219
return err;
11671220
}
@@ -1447,21 +1500,63 @@ static int decoderEngine(decodeRec *dec) {
14471500
err = decodeGrid(dec, i, 0);
14481501

14491502
// *** make sure decode fits somewhere ***
1450-
if (isRestricted(i)) {
1503+
if (err==0 && isRestricted(i)) {
14511504
int fitssomewhere = 0;
14521505
int j;
14531506
for (j = i - 1; j >= from; j--) { // look in previous rects
14541507
if (!isRestricted(j)) {
1455-
if (fitsInsideWithRoom(&dec->coord32, j)) {
1508+
if (fitsInside(&dec->coord32, j)) {
14561509
fitssomewhere = 1;
14571510
break;
14581511
}
14591512
}
14601513
}
1514+
#ifdef FORCE_RECODE
1515+
if (err==0 && !fitssomewhere) {
1516+
for (j = from; j < i; j++) { // try all smaller rectangles j
1517+
if (!isRestricted(j)) {
1518+
const mminforec *b = boundaries(j);
1519+
int bminx = b->minx;
1520+
int bmaxx = b->maxx;
1521+
point p = dec->result;
1522+
1523+
if (bmaxx < 0 && p.lon > 0) {
1524+
bminx += 360000000;
1525+
bmaxx += 360000000;
1526+
} //
1527+
1528+
if (p.lat < b->miny && b->miny <= dec->range_max.lat) {
1529+
p.lat = b->miny;
1530+
} // keep in range
1531+
if (p.lat >= b->maxy && b->maxy > dec->range_min.lat) {
1532+
p.lat = b->maxy - MICROMETER;
1533+
} // keep in range
1534+
if (p.lon < bminx && bminx <= dec->range_max.lon) {
1535+
p.lon = bminx;
1536+
} // keep in range
1537+
if (p.lon >= bmaxx && bmaxx > dec->range_min.lon) {
1538+
p.lon = bmaxx - MICROMETER;
1539+
} // keep in range
1540+
1541+
if ( p.lat > dec->range_min.lat && p.lat < dec->range_max.lat &&
1542+
p.lon > dec->range_min.lon && p.lon < dec->range_max.lon &&
1543+
b->miny <= p.lat && p.lat < b->maxy &&
1544+
bminx <= p.lon && p.lon < bmaxx ) {
1545+
1546+
dec->result = p;
1547+
dec->coord32.lat = (int) floor(p.lat);
1548+
dec->coord32.lon = (int) floor(p.lon);
1549+
fitssomewhere = 1;
1550+
break;
1551+
} // candidate
1552+
} // isRestricted
1553+
} // for j
1554+
}
1555+
#endif //FORCE_RECODE
14611556
if (!fitssomewhere) {
14621557
err = -1234;
14631558
}
1464-
} // *** make sure decode fits somewhere ***
1559+
} // *** make sure decode fits somewhere ***
14651560
break;
14661561
}
14671562
}
@@ -1507,8 +1602,36 @@ static int decoderEngine(decodeRec *dec) {
15071602
if (err == 0) {
15081603
if (ccode != ccode_earth) {
15091604
if (!fitsInsideWithRoom(&dec->coord32, lastrec(ccode))) {
1510-
err = -2222; // EVER?
1605+
err = -2222;
15111606
}
1607+
#ifdef FORCE_RECODE
1608+
else {
1609+
const mminforec *b = boundaries(lastrec(ccode));
1610+
int bminx = b->minx;
1611+
int bmaxx = b->maxx;
1612+
if (dec->coord32.lon < 0 && bminx > 0) {
1613+
bminx -= 360000000;
1614+
bmaxx -= 360000000;
1615+
}
1616+
1617+
if (dec->result.lat < b->miny / 1000000.0) {
1618+
dec->result.lat = b->miny / 1000000.0;
1619+
} // keep in encompassing territory
1620+
if (dec->result.lat >= b->maxy / 1000000.0) {
1621+
dec->result.lat = (b->maxy - MICROMETER) / 1000000.0;
1622+
} // keep in encompassing territory
1623+
1624+
if (dec->result.lon < bminx / 1000000.0) {
1625+
dec->result.lon = bminx / 1000000.0;
1626+
} // keep in encompassing territory
1627+
if (dec->result.lon >= bmaxx / 1000000.0) {
1628+
dec->result.lon = (bmaxx - MICROMETER) / 1000000.0;
1629+
} // keep in encompassing territory
1630+
1631+
dec->coord32.lat = (int) floor(dec->result.lat * 1000000);
1632+
dec->coord32.lon = (int) floor(dec->result.lon * 1000000);
1633+
} // FORCE_RECODE
1634+
#endif
15121635
}
15131636
}
15141637

@@ -1764,37 +1887,36 @@ static int encodeLatLonToMapcodes_internal(char **v, Mapcodes *mapcodes, double
17641887
enc.mapcodes->count = 0;
17651888

17661889
if (lat < -90) { lat = -90; } else if (lat > 90) { lat = 90; }
1767-
if (lon < -180 || lon > 180) {
1768-
lon -= (360.0 * floor(lon / 360));
1769-
if (lon >= 180) { lon -= 360; }
1770-
}
1890+
lat += 90; // lat now [0..180]
1891+
lon -= (360.0 * floor(lon / 360)); // lon now in [0..360>
17711892

1772-
lat += 90;
1773-
lon += 180;
17741893
lat *= 1000000;
17751894
lon *= 1000000;
17761895
enc.coord32.lat = (int) lat;
17771896
enc.coord32.lon = (int) lon;
17781897
enc.fraclat = lat - enc.coord32.lat;
17791898
enc.fraclon = lon - enc.coord32.lon;
1780-
// for 8-digit precision, cells are divided into 810,000 by 810,000 minicells.
1781-
enc.fraclat *= 810000;
1782-
if (enc.fraclat < 1) { enc.fraclat = 0; } else {
1783-
if (enc.fraclat > 809999) {
1784-
enc.fraclat = 0;
1785-
enc.coord32.lat++;
1786-
} else { enc.fraclat /= 810000; }
1787-
}
1788-
enc.fraclon *= 810000;
1789-
if (enc.fraclon < 1) { enc.fraclon = 0; } else {
1790-
if (enc.fraclon > 809999) {
1791-
enc.fraclon = 0;
1792-
enc.coord32.lon++;
1793-
} else { enc.fraclon /= 810000; }
1899+
{
1900+
double f;
1901+
// for 8-digit precision, cells are divided into 810,000 by 810,000 minicells.
1902+
f = enc.fraclat * MAX_PRECISION_FACTOR;
1903+
if (f < 1) { enc.fraclat = 0; } else {
1904+
if (f >= (MAX_PRECISION_FACTOR - 0.5)) {
1905+
enc.fraclat = 0;
1906+
enc.coord32.lat++;
1907+
}
1908+
}
1909+
f = enc.fraclon * MAX_PRECISION_FACTOR;
1910+
if (f < 1) { enc.fraclon = 0; } else {
1911+
if (f >= (MAX_PRECISION_FACTOR - 0.5)) {
1912+
enc.fraclon = 0;
1913+
enc.coord32.lon++;
1914+
}
1915+
}
17941916
}
17951917
enc.coord32.lat -= 90000000;
1796-
enc.coord32.lon %= 360000000;
1797-
enc.coord32.lon -= 180000000;
1918+
if (enc.coord32.lon >= 180000000)
1919+
enc.coord32.lon -= 360000000;
17981920

17991921
if (tc <= 0) // ALL results?
18001922
{

0 commit comments

Comments
 (0)