Patch over bad behaviour for near-parallel boxes and distancetree.

Unfortunately, this "fix" is probably breakable, just put the boxes closer and closer together until we get "close enough to be the same" results from the edge normal comparison again. Only way around this is to do dot product and comparisons in higher-than-double space
References #4223


git-svn-id: http://svn.osgeo.org/postgis/trunk@16981 b70326c6-7e19-0410-871a-916f4a2858ee

Author: pramsey

liblwgeom/cunit/cu_geodetic.c

@@ -586,6 +586,12 @@ static void test_edge_intersects(void)
 	GEOGRAPHIC_POINT g;
 	uint32_t rv;
 
+	/* 5m close case */
+	line2pts("LINESTRING(58.5112113206308 0, 58.511211320077201 0.00090193752520337797)", &A1, &A2);
+	line2pts("LINESTRING(58.511166525601702 0.00027058124084120699, 58.511166525562899 0.00036077498778824899)", &B1, &B2);
+	rv = edge_intersects(&A1, &A2, &B1, &B2);
+	CU_ASSERT(rv == 0);
+
 	/* Covers case, end-to-end intersection */
 	line2pts("LINESTRING(50 -10.999999999999998224, -10.0 50.0)", &A1, &A2);
 	line2pts("LINESTRING(-10.0 50.0, -10.272779983831613393 -16.937003313332997578)", &B1, &B2);

liblwgeom/cunit/cu_tree.c

@@ -177,11 +177,35 @@ static void test_tree_circ_distance(void)
 	LWGEOM *lwg1, *lwg2;
 	CIRC_NODE *c1, *c2;
 	SPHEROID s;
-	double d1, d2, d3, d4;
+	double d1, d2, d3, d4, e1, e2;
 	double threshold = 0.0;
 
 	spheroid_init(&s, 1.0, 1.0);
 
+	/* Ticket #4223 */
+	/* tall skinny rectangle */
+	lwg1 = lwgeom_from_wkt("POLYGON((58.5112113206308 0,58.5112113200772 0.000901937525203378,58.511300910044 0.000901937636668872,58.5113009105976 0,58.5112113206308 0))", LW_PARSER_CHECK_NONE);
+	/* square box 5m to left */
+	lwg2 = lwgeom_from_wkt("POLYGON((58.5111665256017 0.000270581240841207,58.5111665255629 0.000360774987788249,58.5110769356128 0.000360774943200728,58.5110769356515 0.000270581207400566,58.5111665256017 0.000270581240841207))", LW_PARSER_CHECK_NONE);
+	c1 = lwgeom_calculate_circ_tree(lwg1);
+	c2 = lwgeom_calculate_circ_tree(lwg2);
+	d1 = circ_tree_distance_tree(c1, c2, &s, threshold);
+	d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold);
+	e1 = d1 * WGS84_RADIUS;
+	e2 = d2 * WGS84_RADIUS;
+	// printf("d1 = %g   d2 = %g\n", d1, d2);
+	// printf("e1 = %g   e2 = %g\n", e1, e2);
+	// printf("polygon a\n");
+	// circ_tree_print(c1, 0);
+	// printf("polygon b\n");
+	// circ_tree_print(c2, 0);
+	circ_tree_free(c1);
+	circ_tree_free(c2);
+	lwgeom_free(lwg1);
+	lwgeom_free(lwg2);
+	CU_ASSERT_DOUBLE_EQUAL(e1, e2, 0.0001);
+
+
 	/* Ticket #1958 */
 	lwg1 = lwgeom_from_wkt("LINESTRING(22.88333 41.96667,21.32667 42.13667)", LW_PARSER_CHECK_NONE);
 	lwg2 = lwgeom_from_wkt("POLYGON((22.94472 41.34667,22.87528 41.99028,22.87389 41.98472,22.87472 41.98333,22.94472 41.34667))", LW_PARSER_CHECK_NONE);

liblwgeom/lwgeodetic.h

@@ -37,6 +37,11 @@
 #define NAN 0.0/0.0
 #endif
 
+#ifdef FP_TOLERANCE
+#undef FP_TOLERANCE
+#define FP_TOLERANCE 1e-14
+#endif 
+
 extern int gbox_geocentric_slow;
 
 #define POW2(x) ((x)*(x))