#2018, Distance calculation support for arc features (circstring, com…

…poundcurve, curvepolygon)

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

NEWS

@@ -59,6 +59,8 @@ PostGIS 2.1.0
   - Casts to/from PostgreSQL geotypes (point/path/polygon).
   - #310, ST_DumpPoints converted to a C function (Nathan Wagner)
   - #2011, ST_DumpValues to output raster as array (Bborie Park / UC Davis)
+  - #2018, ST_Distance support for CircularString, CurvePolygon, MultiCurve,
+    MultiSurface, CompoundCurve 
   - #2030, n-raster (and n-band) ST_MapAlgebra (Bborie Park / UC Davis)
   - Added geomval array variant of ST_SetValues() to set many pixel values of
     a band using a set of geometries and corresponding values in one call

liblwgeom/cunit/cu_measures.c

@@ -28,82 +28,167 @@ static LWGEOM* lwgeom_from_text(const char *str)
 	return r.geom;
 }
 
-static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res)
+#define DIST2DTEST(str1, str2, res) do_test_mindistance2d_tolerance(str1, str2, res, __LINE__)
+
+static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line)
 {
 	LWGEOM *lw1;
 	LWGEOM *lw2;
 	double distance;
+	char *msg1 = "test_mindistance2d_tolerance failed (got %g expected %g) at line %d\n";
+	char *msg2 = "\n\ndo_test_mindistance2d_tolerance: NULL lwgeom generated from WKT\n  %s\n\n";
 
 	lw1 = lwgeom_from_wkt(in1, LW_PARSER_CHECK_NONE);
 	lw2 = lwgeom_from_wkt(in2, LW_PARSER_CHECK_NONE);
 
+	if ( ! lw1 )
+	{
+		printf(msg2, in1);
+		exit(1);
+	}
+	if ( ! lw2 )
+	{
+		printf(msg2, in2);
+		exit(1);
+	}
+
 	distance = lwgeom_mindistance2d_tolerance(lw1, lw2, 0.0);
-	CU_ASSERT_EQUAL(distance, expected_res);
 	lwgeom_free(lw1);
 	lwgeom_free(lw2);
 
+	if ( fabs(distance - expected_res) > 0.00001 )
+	{
+		printf(msg1, distance, expected_res, line);
+		CU_FAIL();
+	}
+	else
+	{
+		CU_PASS();
+	}
+
 }
+
 static void test_mindistance2d_tolerance(void)
 {
 	/*
 	** Simple case.
 	*/
-	do_test_mindistance2d_tolerance("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
+	DIST2DTEST("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
 
 	/*
 	** Point vs Geometry Collection.
 	*/
-	do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
+	DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
 
 	/*
 	** Point vs Geometry Collection Collection.
 	*/
-	do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
+	DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
 
 	/*
 	** Point vs Geometry Collection Collection Collection.
 	*/
-	do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
+	DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
 
 	/*
 	** Point vs Geometry Collection Collection Collection Multipoint.
 	*/
-	do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
+	DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
 
 	/*
 	** Geometry Collection vs Geometry Collection
 	*/
-	do_test_mindistance2d_tolerance("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
+	DIST2DTEST("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
 
 	/*
 	** Geometry Collection Collection vs Geometry Collection Collection
 	*/
-	do_test_mindistance2d_tolerance("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
+	DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
 
 	/*
 	** Geometry Collection Collection Multipoint vs Geometry Collection Collection Multipoint
 	*/
-	do_test_mindistance2d_tolerance("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0);
+	DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0);
 
 	/*
 	** Linestring vs its start point 
 	*/
-	do_test_mindistance2d_tolerance("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
+	DIST2DTEST("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
 
 	/*
 	** Linestring vs its end point 
 	*/
-	do_test_mindistance2d_tolerance("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
+	DIST2DTEST("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
 
 	/*
 	** Linestring vs its start point (tricky number, see #1459)
 	*/
-	do_test_mindistance2d_tolerance("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
+	DIST2DTEST("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
 
 	/*
 	** Linestring vs its end point (tricky number, see #1459)
 	*/
-	do_test_mindistance2d_tolerance("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
+	DIST2DTEST("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
+
+	/*
+	* Circular string and point 
+	*/
+	DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "POINT(0 0)", 1);
+	DIST2DTEST("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)", "POINT(0 0)", 1);
+
+	/*
+	* Circular string and Circular string 
+	*/
+	DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "CIRCULARSTRING(0 0, 1 -1, 2 0)", 1);
+
+	/*
+	* CurvePolygon and Point
+	*/
+	static char *cs1 = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))";
+	DIST2DTEST(cs1, "POINT(3 14)", 1);
+	DIST2DTEST(cs1, "POINT(3 8)", 0);
+	DIST2DTEST(cs1, "POINT(6 5)", 1);
+	DIST2DTEST(cs1, "POINT(6 4)", 0);
+
+	/*
+	* CurvePolygon and Linestring
+	*/
+	DIST2DTEST(cs1, "LINESTRING(0 0, 50 0)", 0.917484);
+	DIST2DTEST(cs1, "LINESTRING(6 0, 10 7)", 0);
+	DIST2DTEST(cs1, "LINESTRING(4 4, 4 8)", 0);
+	DIST2DTEST(cs1, "LINESTRING(4 7, 5 6, 6 7)", 0.585786);
+	DIST2DTEST(cs1, "LINESTRING(10 0, 10 2, 10 0)", 1.52913);
+
+	/*
+	* CurvePolygon and Polygon
+	*/
+	DIST2DTEST(cs1, "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))", 0.58415);
+	DIST2DTEST(cs1, "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))", 0);
+	DIST2DTEST(cs1, "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))", 0);
+
+	/*
+	* CurvePolygon and CurvePolygon
+	*/
+	DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))", 0.0475666);
+	DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))", 0.0);
+
+	/*
+	* MultiSurface and CurvePolygon 
+	*/
+	static char *cs2 = "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))";
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
+	DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
+
+	/*
+	* MultiCurve and Linestring
+	*/
+	DIST2DTEST("LINESTRING(0.5 1,0.5 3)", "MULTICURVE(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3),(0 0, 0 5))", 0.5);
 
 }
 

liblwgeom/cunit/cu_ptarray.c

@@ -409,7 +409,7 @@ static void test_ptarray_desegmentize()
 
 static void test_ptarray_contains_point() 
 {
-/* int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) */
+/* int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt, int *winding_number) */
 
 	LWLINE *lwline;
 	POINTARRAY *pa;
@@ -496,119 +496,140 @@ static void test_ptarray_contains_point()
 	lwline_free(lwline);
 }
 
-
-static void test_ptarray_contains_point_arc() 
+static void test_ptarrayarc_contains_point() 
 {
-	/* int ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) */
+	/* int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt) */
 
 	LWLINE *lwline;
 	POINTARRAY *pa;
 	POINT2D pt;
 	int rv;
 
-	/* Collection of semi-circles surrounding unit square */
+	/*** Collection of semi-circles surrounding unit square ***/
 	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
 	pa = lwline->points;
 
 	/* Point in middle of square */
 	pt.x = 0;
 	pt.y = 0;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_INSIDE);
 
 	/* Point in left lobe */
 	pt.x = -1.1;
 	pt.y = 0.1;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_INSIDE);	
 
 	/* Point on boundary of left lobe */
 	pt.x = -1;
 	pt.y = 0;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_INSIDE);	
 
 	/* Point on boundary vertex */
 	pt.x = -1;
 	pt.y = 1;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_BOUNDARY);	
 
 	/* Point outside */
 	pt.x = -1.5;
 	pt.y = 1.5;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_OUTSIDE);	
 
-	/* Two-edge ring made up of semi-circles (really, a circle) */
+	/*** Two-edge ring made up of semi-circles (really, a circle) ***/
 	lwline_free(lwline);
-
 	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0, 0 -1, -1 0)"));
 	pa = lwline->points;
 
 	/* Point outside */
 	pt.x = -1.5;
 	pt.y = 1.5;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
+	CU_ASSERT_EQUAL(rv, LW_OUTSIDE);	
+
+	/* Point more outside */
+	pt.x = 2.5;
+	pt.y = 1.5;
+	rv = ptarrayarc_contains_point(pa, &pt);
+	CU_ASSERT_EQUAL(rv, LW_OUTSIDE);	
+
+	/* Point more outside */
+	pt.x = 2.5;
+	pt.y = 2.5;
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_OUTSIDE);	
 
 	/* Point inside at middle */
 	pt.x = 0;
 	pt.y = 0;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_INSIDE);	
 
 	/* Point inside offset from middle */
 	pt.x = 0.01;
 	pt.y = 0.01;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_INSIDE);	
 
 	/* Point on edge vertex */
 	pt.x = 0;
 	pt.y = 1;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_BOUNDARY);	
 
-	/* One-edge ring, closed circle */
+	/*** Two-edge ring, closed ***/
+	lwline_free(lwline);
+	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1 6, 6 1, 9 7, 6 10, 1 6)"));
+	pa = lwline->points;
+
+	/* Point to left of ring */
+	pt.x = 20;
+	pt.y = 4;
+	rv = ptarrayarc_contains_point(pa, &pt);
+	CU_ASSERT_EQUAL(rv, LW_OUTSIDE);	
+
+	/*** One-edge ring, closed circle ***/
 	lwline_free(lwline);
 	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, -1 0)"));
 	pa = lwline->points;
 
 	/* Point inside */
 	pt.x = 0;
 	pt.y = 0;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_INSIDE);	
 
 	/* Point outside */
 	pt.x = 0;
 	pt.y = 2;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_OUTSIDE);	
 
 	/* Point on boundary */
 	pt.x = 0;
 	pt.y = 1;
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	CU_ASSERT_EQUAL(rv, LW_BOUNDARY);	
 
-	/* Overshort ring */
+	/*** Overshort ring ***/
 	lwline_free(lwline);
 	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0)"));
 	pa = lwline->points;
 	cu_error_msg_reset();
-	rv = ptarray_contains_point_arc(pa, &pt);
+	rv = ptarrayarc_contains_point(pa, &pt);
 	//printf("%s\n", cu_error_msg);
-	CU_ASSERT_STRING_EQUAL("ptarray_contains_point_arc called with even number of points", cu_error_msg);
+	CU_ASSERT_STRING_EQUAL("ptarrayarc_contains_point called with even number of points", cu_error_msg);
 
-	/* Unclosed ring */
+	/*** Unclosed ring ***/
 	lwline_free(lwline);
-	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, 1 0)"));
+	lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, 2 0)"));
 	pa = lwline->points;
 	cu_error_msg_reset();
-	rv = ptarray_contains_point_arc(pa, &pt);
-	CU_ASSERT_STRING_EQUAL("ptarray_contains_point_arc called on unclosed ring", cu_error_msg);
+	rv = ptarrayarc_contains_point(pa, &pt);
+	CU_ASSERT_STRING_EQUAL("ptarrayarc_contains_point called on unclosed ring", cu_error_msg);
 
 	lwline_free(lwline);
 }
@@ -627,7 +648,7 @@ CU_TestInfo ptarray_tests[] =
 	PG_TEST(test_ptarray_desegmentize),
 	PG_TEST(test_ptarray_insert_point),
 	PG_TEST(test_ptarray_contains_point),
-	PG_TEST(test_ptarray_contains_point_arc),
+	PG_TEST(test_ptarrayarc_contains_point),
 	CU_TEST_INFO_NULL
 };
 CU_SuiteInfo ptarray_suite = {"ptarray", NULL, NULL, ptarray_tests };

liblwgeom/liblwgeom.h.in

@@ -919,9 +919,11 @@ extern int lwcurvepoly_add_ring(LWCURVEPOLY *poly, LWGEOM *ring);
 */
 extern int lwcompound_add_lwgeom(LWCOMPOUND *comp, LWGEOM *geom);
 
-
-
-
+/**
+* Construct an equivalent curve polygon from a polygon. Curve polygons
+* can have linear rings as their rings, so this works fine (in theory?)
+*/
+extern LWCURVEPOLY* lwcurvepoly_construct_from_lwpoly(LWPOLY *lwpoly);
 
 
 /******************************************************************