Improve precision of polygon area calculations (#19)

* Improve precision of get_area by using long double

* Trying to get consistent polygon area results between ARM and x86

* Calculate polygon area closer to the origin for better precision

* Update changelog

* Also exercise tiny polygon dust in the ring area test

They previously behaved differently here between x86 and ARM

* On M1 Macs, long double is just double anyway, so don't use it

* Be more careful about overflow: scale the polygon ring down into range

* Fix the bug I just introduced in the scaled area calculation

* Use only the sign from the scaled-down area calculation

Co-authored-by: Roman Karavia <[email protected]>

CHANGELOG.md

@@ -1,3 +1,7 @@
+## 2.8.1
+
+* Improve precision of polygon ring area calculations
+
 ## 2.8.0
 
 * Add the option to use a different simplification level at maxzoom

decode.cpp

@@ -207,9 +207,9 @@ void handle(std::string message, int z, unsigned x, unsigned y, std::set<std::st
 		}
 
 		double scale = 0;
-		if (coordinate_mode == 1) { // fraction
+		if (coordinate_mode == 1) {  // fraction
 			scale = layer.extent;
-		} else if (coordinate_mode == 2) { // integer
+		} else if (coordinate_mode == 2) {  // integer
 			scale = 1;
 		}
 		layer_to_geojson(layer, z, x, y, !pipeline, pipeline, pipeline, false, 0, 0, 0, !force, state, scale);

geometry.cpp

@@ -160,13 +160,88 @@ drawvec remove_noop(drawvec geom, int type, int shift) {
 	return out;
 }
 
-double get_area(drawvec &geom, size_t i, size_t j) {
+double get_area_scaled(const drawvec &geom, size_t i, size_t j) {
+	const double max_exact_double = (double) ((1LL << 53) - 1);
+
+	// keep scaling the geometry down until we can calculate its area without overflow
+	for (long long scale = 2; scale < (1LL << 30); scale *= 2) {
+		long long bx = geom[i].x;
+		long long by = geom[i].y;
+		bool again = false;
+
+		// https://en.wikipedia.org/wiki/Shoelace_formula
+		double area = 0;
+		for (size_t k = i; k < j; k++) {
+			area += (double) ((geom[k].x - bx) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].y - by) / scale);
+			if (std::fabs(area) >= max_exact_double) {
+				again = true;
+				break;
+			}
+			area -= (double) ((geom[k].y - by) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].x - bx) / scale);
+			if (std::fabs(area) >= max_exact_double) {
+				again = true;
+				break;
+			}
+		}
+
+		if (again) {
+			continue;
+		} else {
+			area /= 2;
+			return area * scale * scale;
+		}
+	}
+
+	fprintf(stderr, "get_area_scaled: can't happen\n");
+	exit(EXIT_IMPOSSIBLE);
+}
+
+double get_area(const drawvec &geom, size_t i, size_t j) {
+	const double max_exact_double = (double) ((1LL << 53) - 1);
+
+	// Coordinates in `geom` are 40-bit integers, so there is no good way
+	// to multiply them without possible precision loss. Since they probably
+	// do not use the full precision, shift them nearer to the origin so
+	// their product is more likely to be exactly representable as a double.
+	//
+	// (In practice they are actually 34-bit integers: 32 bits for the
+	// Mercator world plane, plus another two bits so features can stick
+	// off either the left or right side. But that is still too many bits
+	// for the product to fit either in a 64-bit long long or in a
+	// double where the largest exact integer is 2^53.)
+	//
+	// If the intermediate calculation still exceeds 2^53, start trying to
+	// recalculate the area by scaling down the geometry. This will not
+	// produce as precise an area, but it will still be close, and the
+	// sign will be correct, which is more important, since the sign
+	// determines the winding order of the rings. We can then use that
+	// sign with this generally more precise area calculation.
+
+	long long bx = geom[i].x;
+	long long by = geom[i].y;
+
+	// https://en.wikipedia.org/wiki/Shoelace_formula
 	double area = 0;
+	bool overflow = false;
 	for (size_t k = i; k < j; k++) {
-		area += (long double) geom[k].x * (long double) geom[i + ((k - i + 1) % (j - i))].y;
-		area -= (long double) geom[k].y * (long double) geom[i + ((k - i + 1) % (j - i))].x;
+		area += (double) (geom[k].x - bx) * (double) (geom[i + ((k - i + 1) % (j - i))].y - by);
+		if (std::fabs(area) >= max_exact_double) {
+			overflow = true;
+		}
+		area -= (double) (geom[k].y - by) * (double) (geom[i + ((k - i + 1) % (j - i))].x - bx);
+		if (std::fabs(area) >= max_exact_double) {
+			overflow = true;
+		}
 	}
 	area /= 2;
+
+	if (overflow) {
+		double scaled_area = get_area_scaled(geom, i, j);
+		if ((area < 0 && scaled_area > 0) || (area > 0 && scaled_area < 0)) {
+			area = -area;
+		}
+	}
+
 	return area;
 }
 
@@ -518,7 +593,7 @@ drawvec simple_clip_poly(drawvec &geom, int z, int buffer) {
 
 drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double *accum_area) {
 	drawvec out;
-	const long long pixel = (1 << (32 - detail - z)) * tiny_polygon_size;
+	const double pixel = (1LL << (32 - detail - z)) * (double) tiny_polygon_size;
 
 	*reduced = true;
 	bool included_last_outer = false;

geometry.hpp

@@ -74,7 +74,7 @@ drawvec reorder_lines(drawvec &geom);
 drawvec fix_polygon(drawvec &geom);
 std::vector<drawvec> chop_polygon(std::vector<drawvec> &geoms);
 void check_polygon(drawvec &geom);
-double get_area(drawvec &geom, size_t i, size_t j);
+double get_area(const drawvec &geom, size_t i, size_t j);
 double get_mp_area(drawvec &geom);
 
 drawvec simple_clip_poly(drawvec &geom, long long x1, long long y1, long long x2, long long y2);

main.cpp

@@ -2289,7 +2289,7 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 				if (maxzoom == 0) {
 					droprate = 2.5;
 				} else {
-					droprate = exp(log((long double) max[0].count / max[maxzoom].count) / (maxzoom));
+					droprate = exp(log((double) max[0].count / max[maxzoom].count) / (maxzoom));
 					if (!quiet) {
 						fprintf(stderr, "Choosing a drop rate of -r%f to get from %lld to %lld in %d zooms\n", droprate, max[maxzoom].count, max[0].count, maxzoom);
 					}
@@ -2298,7 +2298,7 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 
 			basezoom = 0;
 			for (z = 0; z <= maxzoom; z++) {
-				double zoomdiff = log((long double) max[z].count / max_features) / log(droprate);
+				double zoomdiff = log((double) max[z].count / max_features) / log(droprate);
 				if (zoomdiff + z > basezoom) {
 					basezoom = ceil(zoomdiff + z);
 				}
@@ -2314,7 +2314,7 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 				double interval = exp(log(droprate) * (basezoom - z));
 
 				if (max[z].count / interval >= max_features) {
-					interval = (long double) max[z].count / max_features;
+					interval = (double) max[z].count / max_features;
 					droprate = exp(log(interval) / (basezoom - z));
 					interval = exp(log(droprate) * (basezoom - z));
 

tests/single-polygons/in.json

@@ -0,0 +1,30 @@
+{"type":"Feature","properties":{"seq":6},"geometry":{"type":"Polygon","coordinates":[[[174.72151704933612,-37.16620984485238],[174.721517058737,-37.16620953185807],[174.7215172151761,-37.16620953486748],[174.72151720953557,-37.16620972266407],[174.72151705121628,-37.166209782253527],[174.72151704933612,-37.16620984485238]]]}}
+{"type":"Feature","properties":{"seq":7},"geometry":{"type":"Polygon","coordinates":[[[174.7215178334118,-37.16620979730054],[174.72151760063333,-37.16620973018757],[174.72151760627387,-37.16620954239099],[174.72151784093254,-37.16620954690509],[174.72151791351156,-37.166209736206379],[174.7215178334118,-37.16620979730054]]]}}
+{"type":"Feature","properties":{"seq":8},"geometry":{"type":"Polygon","coordinates":[[[174.72151861936775,-37.16620968714982],[174.7215183064895,-37.166209681131018],[174.72151823203033,-37.16620955442859],[174.7215186231281,-37.1662095619521],[174.72151861936775,-37.16620968714982]]]}}
+{"type":"Feature","properties":{"seq":9},"geometry":{"type":"Polygon","coordinates":[[[174.72151955424205,-37.16620983040395],[174.72151932710407,-37.1662095754944],[174.7215196399823,-37.16620958151319],[174.72151963434178,-37.16620976930978],[174.72151955424205,-37.16620983040395]]]}}
+{"type":"Feature","properties":{"seq":10},"geometry":{"type":"Polygon","coordinates":[[[174.7215202600982,-37.166209781347379],[174.7215201036591,-37.16620977833798],[174.72151995286053,-37.166209587532],[174.7215203439583,-37.16620959505549],[174.72152034019795,-37.166209720253217],[174.7215202600982,-37.166209781347379]]]}}
+{"type":"Feature","properties":{"seq":11},"geometry":{"type":"Polygon","coordinates":[[[174.72152111863316,-37.16620986049794],[174.7215208076351,-37.16620979188027],[174.72152081327563,-37.16620960408368],[174.72152112615385,-37.16620961010249],[174.72152111863316,-37.16620986049794]]]}}
+{"type":"Feature","properties":{"seq":12},"geometry":{"type":"Polygon","coordinates":[[[174.72152182448932,-37.16620981144135],[174.72152151349128,-37.166209742823699],[174.7215215172516,-37.166209617625977],[174.7215219083494,-37.16620962514947],[174.72152190458906,-37.16620975034719],[174.72152182448932,-37.16620981144135]]]}}
+{"type":"Feature","properties":{"seq":13},"geometry":{"type":"Polygon","coordinates":[[[174.7215236235391,-37.1662098460494],[174.72152331630142,-37.16620965223402],[174.72152370739918,-37.166209659757509],[174.7215237036388,-37.16620978495524],[174.7215236235391,-37.1662098460494]]]}}
+{"type":"Feature","properties":{"seq":14},"geometry":{"type":"Polygon","coordinates":[[[174.72152440761483,-37.16620979849751],[174.72152409849697,-37.16620966728101],[174.7215244895947,-37.16620967480448],[174.721524409495,-37.166209735898657],[174.72152440761483,-37.16620979849751]]]}}
+{"type":"Feature","properties":{"seq":15},"geometry":{"type":"Polygon","coordinates":[[[174.72152526802993,-37.16620981504918],[174.72152495703188,-37.166209746431537],[174.72152527179029,-37.16620968985145],[174.72152526802993,-37.16620981504918]]]}}
+{"type":"Feature","properties":{"seq":16},"geometry":{"type":"Polygon","coordinates":[[[174.72152581180647,-37.16620995077977],[174.72152565536738,-37.16620994777038],[174.72152558466849,-37.166209695870239],[174.7215258975467,-37.16620970188902],[174.7215258919062,-37.1662098896856],[174.72152581180647,-37.16620995077977]]]}}
+{"type":"Feature","properties":{"seq":17},"geometry":{"type":"Polygon","coordinates":[[[174.72152690876045,-37.16620990924664],[174.72152659776237,-37.166209840629],[174.72152660152273,-37.16620971543128],[174.72152691440093,-37.16620972145006],[174.72152690876045,-37.16620990924664]]]}}
+{"type":"Feature","properties":{"seq":18},"geometry":{"type":"Polygon","coordinates":[[[174.72153363564224,-37.16621003865021],[174.72153324642464,-37.1662099685279],[174.72153325018497,-37.166209843330168],[174.72153364128273,-37.16620985085363],[174.72153363564224,-37.16621003865021]]]}}
+{"type":"Feature","properties":{"seq":19},"geometry":{"type":"Polygon","coordinates":[[[174.72153504359427,-37.166210065734649],[174.72153488715515,-37.16621006272526],[174.72153473635653,-37.166209871919289],[174.72153520567387,-37.166209880947437],[174.72153512369398,-37.16621000464046],[174.72153504359427,-37.166210065734649]]]}}
+{"type":"Feature","properties":{"seq":20},"geometry":{"type":"Polygon","coordinates":[[[174.7215366079854,-37.166210095828429],[174.72153621876778,-37.16621002570611],[174.7215362225281,-37.16620990050839],[174.72153661362587,-37.166209908031827],[174.7215366079854,-37.166210095828429]]]}}
+{"type":"Feature","properties":{"seq":21},"geometry":{"type":"Polygon","coordinates":[[[174.72153817425667,-37.16621006332332],[174.72153786137847,-37.16621005730457],[174.7215377869192,-37.16620993060216],[174.721538178017,-37.16620993812559],[174.72153817425667,-37.16621006332332]]]}}
+{"type":"Feature","properties":{"seq":22},"geometry":{"type":"Polygon","coordinates":[[[174.72154004964589,-37.166210162034669],[174.7215398149872,-37.1662101575206],[174.7215397424081,-37.16620996821933],[174.7215401335059,-37.16620997574277],[174.72154012974557,-37.16621010094049],[174.72154004964589,-37.166210162034669]]]}}
+{"type":"Feature","properties":{"seq":23},"geometry":{"type":"Polygon","coordinates":[[[174.72154153581747,-37.166210190623697],[174.7215412248194,-37.16621012200609],[174.7215412285797,-37.16620999680837],[174.72154154145793,-37.16621000282711],[174.72154153581747,-37.166210190623697]]]}}
+{"type":"Feature","properties":{"seq":24},"geometry":{"type":"Polygon","coordinates":[[[174.72154231801304,-37.16621020567055],[174.7215421615739,-37.16621020266118],[174.72154201077528,-37.16621001185522],[174.7215424800926,-37.16621002088333],[174.72154239811273,-37.16621014457637],[174.72154231801304,-37.16621020567055]]]}}
+{"type":"Feature","properties":{"seq":25},"geometry":{"type":"Polygon","coordinates":[[[174.7215440332028,-37.1662104265702],[174.7215438041846,-37.16621023425956],[174.72154349130637,-37.16621022824081],[174.72154341872727,-37.16621003893954],[174.72154380982509,-37.16621004646295],[174.72154380606475,-37.16621017166068],[174.72154404072342,-37.16621017617474],[174.7215441133025,-37.166210365476008],[174.7215440332028,-37.1662104265702]]]}}
+{"type":"Feature","properties":{"seq":26},"geometry":{"type":"Polygon","coordinates":[[[174.72154474093913,-37.16621031491462],[174.7215445845,-37.16621031190526],[174.7215445138011,-37.166210060005109],[174.7215448266793,-37.166210066023847],[174.72154482103884,-37.16621025382044],[174.72154474093913,-37.16621031491462]]]}}
+{"type":"Feature","properties":{"seq":27},"geometry":{"type":"Polygon","coordinates":[[[174.7215463760292,-37.1662105969084],[174.72154606503114,-37.16621052829082],[174.7215461564118,-37.16621009160346],[174.72154639107047,-37.166210096117499],[174.72154637978952,-37.16621047171068],[174.7215463760292,-37.1662105969084]]]}}
+{"type":"Feature","properties":{"seq":28},"geometry":{"type":"Polygon","coordinates":[[[174.7215470093063,-37.166210358550419],[174.72154677652777,-37.1662102914375],[174.72154678216823,-37.166210103640917],[174.72154709504648,-37.166210109659647],[174.721547089406,-37.16621029745623],[174.7215470093063,-37.166210358550419]]]}}
+{"type":"Feature","properties":{"seq":29},"geometry":{"type":"Polygon","coordinates":[[[174.72154865379714,-37.16621032754985],[174.7215484191385,-37.16621032303581],[174.7215482683398,-37.16621013222985],[174.72154873765718,-37.16621014125794],[174.72154873389688,-37.16621026645566],[174.72154865379714,-37.16621032754985]]]}}
+{"type":"Feature","properties":{"seq":30},"geometry":{"type":"Polygon","coordinates":[[[174.7215495123321,-37.1662104067002],[174.72154920321419,-37.16621027548374],[174.7215492069745,-37.16621015028602],[174.72154951985272,-37.166210156304739],[174.7215495123321,-37.1662104067002]]]}}
+{"type":"Feature","properties":{"seq":31},"geometry":{"type":"Polygon","coordinates":[[[174.7215563193138,-37.166210475008309],[174.7215560846551,-37.166210470494288],[174.72155609029557,-37.16621028269769],[174.7215563249542,-37.166210287211708],[174.7215563193138,-37.166210475008309]]]}}
+{"type":"Feature","properties":{"seq":32},"geometry":{"type":"Polygon","coordinates":[[[174.72155686497056,-37.1662105481399],[174.72155671041157,-37.16621048253168],[174.72155671417188,-37.16621035733396],[174.72155702893026,-37.1662103007538],[174.72155694507024,-37.1662104870457],[174.72155686497056,-37.1662105481399]]]}}
+{"type":"Feature","properties":{"seq":33},"geometry":{"type":"Polygon","coordinates":[[[174.72155858768097,-37.166210518643868],[174.7215583530223,-37.16621051412985],[174.7215582804432,-37.16621032482858],[174.72155867154096,-37.166210332351969],[174.7215585895611,-37.16621045604501],[174.72155858768097,-37.166210518643868]]]}}
+{"type":"Feature","properties":{"seq":34},"geometry":{"type":"Polygon","coordinates":[[[174.72155944057554,-37.166210785590759],[174.72155920591687,-37.16621078107673],[174.7215592153176,-37.16621046808239],[174.72155906263874,-37.16621033987532],[174.721559375517,-37.16621034589402],[174.72155944621597,-37.16621059779416],[174.72155944057554,-37.166210785590759]]]}}
+{"type":"Feature","properties":{"seq":35},"geometry":{"type":"Polygon","coordinates":[[[174.72156406305005,-37.16621062397093],[174.7215637501718,-37.166210617952248],[174.72156375581225,-37.16621043015566],[174.72156406869048,-37.16621043617434],[174.72156406305005,-37.16621062397093]]]}}