More attempts to solve failures to simplify consistently

geometry.cpp

@@ -993,7 +993,7 @@ static long long square_distance_from_line(long long point_x, long long point_y,
 	long long p2x = segB_x - segA_x;
 	long long p2y = segB_y - segA_y;
 	long long something = p2x * p2x + p2y * p2y;
-	long long u = 0 == something ? 0 : ((point_x - segA_x) * p2x + (point_y - segA_y) * p2y) * scale / something;
+	double u = 0 == something ? 0 : ((point_x - segA_x) * p2x + (point_y - segA_y) * p2y) * scale / something;
 
 	if (u > scale) {
 		u = scale;
@@ -1014,19 +1014,20 @@ static long long square_distance_from_line(long long point_x, long long point_y,
 static void douglas_peucker(drawvec &geom, int start, int n, double e, size_t kept, size_t retain) {
 	e = e * e;
 	std::stack<int> recursion_stack;
+    printf("doug of %d\n", n);
 
-	{
-		int left_border = 0;
-		int right_border = 1;
-		// Sweep linerarily over array and identify those ranges that need to be checked
-		do {
-			if (geom[start + right_border].necessary) {
-				recursion_stack.push(left_border);
-				recursion_stack.push(right_border);
-				left_border = right_border;
-			}
-			++right_border;
-		} while (right_border < n);
+	if (!geom[start + 0].necessary || !geom[start + n - 1].necessary) {
+		fprintf(stderr, "endpoints not marked necessary\n");
+		exit(EXIT_IMPOSSIBLE);
+	}
+
+	int prev = 0;
+	for (int here = 1; here < n; here++) {
+		if (geom[start + here].necessary) {
+			recursion_stack.push(prev);
+			recursion_stack.push(here);
+            prev = here;
+		}
 	}
 
 	while (!recursion_stack.empty()) {
@@ -1036,6 +1037,8 @@ static void douglas_peucker(drawvec &geom, int start, int n, double e, size_t ke
 		int first = recursion_stack.top();
 		recursion_stack.pop();
 
+        printf("sub-doug of %d\n", second - first + 1);
+
 		long long max_distance = -1;
 		int farthest_element_index;
 		if (geom[start + first] < geom[start + second]) {
@@ -1078,10 +1081,14 @@ static void douglas_peucker(drawvec &geom, int start, int n, double e, size_t ke
 			if (1 < farthest_element_index - first) {
 				recursion_stack.push(first);
 				recursion_stack.push(farthest_element_index);
+
+                printf("split1: %d to %d, %d\n", first, farthest_element_index, farthest_element_index - first + 1);
 			}
 			if (1 < second - farthest_element_index) {
 				recursion_stack.push(farthest_element_index);
 				recursion_stack.push(second);
+
+                printf("split2: %d to %d, %d\n", farthest_element_index, second, second - farthest_element_index + 1);
 			}
 		}
 	}
@@ -1322,7 +1329,7 @@ drawvec fix_polygon(drawvec &geom) {
 				long long xd = ring[a].x - xtotal;
 				long long yd = ring[a].y - ytotal;
 				long long d2 = xd * xd + yd * yd;
-				if (d2 > dist2) {
+				if (d2 > dist2 || (d2 == dist2 && ring[a] < ring[furthest])) {
 					dist2 = d2;
 					furthest = a;
 				}
@@ -1335,7 +1342,7 @@ drawvec fix_polygon(drawvec &geom) {
 				long long xd = ring[a].x - ring[furthest].x;
 				long long yd = ring[a].y - ring[furthest].y;
 				long long d2 = xd * xd + yd * yd;
-				if (d2 > dist2b) {
+				if (d2 > dist2b || (d2 == dist2b && ring[a] < ring[furthestb])) {
 					dist2b = d2;
 					furthestb = a;
 				}

tile.cpp

@@ -2309,16 +2309,17 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 									// would keep next.
 
 									far = 0;
-									double which2 = i;
+									size_t which2 = i;
 
 									for (size_t k = i + 1; k < j - 1; k++) {
 										double xd = sf.geometry[k].x - sf.geometry[i].x;
 										double yd = sf.geometry[k].y - sf.geometry[i].y;
 										double xd2 = sf.geometry[k].x - sf.geometry[which].x;
 										double yd2 = sf.geometry[k].y - sf.geometry[which].y;
-										double d = xd * xd + yd * yd + xd2 * xd2 + yd2 * yd2;
-										if (d > far) {
-											far = d;
+										double d1 = xd * xd + yd * yd;
+										double d2 = xd2 * xd2 + yd2 * yd2;
+										if (d1 + d2 > far && d1 != 0 && d2 != 0) {
+											far = d1 + d2;
 											which2 = k;
 										}
 									}