[postgis-tickets] r16952 - ST_LocateBetween[Elevations]: support POLYGON

Darafei komzpa at gmail.com
Thu Oct 25 04:16:03 PDT 2018


Author: komzpa
Date: 2018-10-25 16:16:02 -0700 (Thu, 25 Oct 2018)
New Revision: 16952

Modified:
   trunk/NEWS
   trunk/liblwgeom/cunit/cu_algorithm.c
   trunk/liblwgeom/liblwgeom_internal.h
   trunk/liblwgeom/lwlinearreferencing.c
Log:
ST_LocateBetween[Elevations]: support POLYGON

Output may be invalid for non-convex, area is preserved.

References #4155
Closes https://github.com/postgis/postgis/pull/319



Modified: trunk/NEWS
===================================================================
--- trunk/NEWS	2018-10-25 14:25:19 UTC (rev 16951)
+++ trunk/NEWS	2018-10-25 23:16:02 UTC (rev 16952)
@@ -27,7 +27,7 @@
     Praliaskouski)
   - #3457, Fix raster envelope shortcut in ST_Clip (Sai-bot)
   - #4215, Use floating point compare in ST_DumpAsPolygons (Darafei Praliaskouski)
-  - #4155, Support for GEOMETRYCOLLECTION in ST_LocateBetween (Darafei
+  - #4155, Support for GEOMETRYCOLLECTION and POLYGON in ST_LocateBetween (Darafei
     Praliaskouski)
 
 PostGIS 2.5.0

Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c	2018-10-25 14:25:19 UTC (rev 16951)
+++ trunk/liblwgeom/cunit/cu_algorithm.c	2018-10-25 23:16:02 UTC (rev 16952)
@@ -546,15 +546,15 @@
 static void test_lwline_clip(void)
 {
 	LWCOLLECTION *c;
-	LWLINE *line = NULL;
-	LWLINE *l51 = NULL;
+	LWGEOM *line = NULL;
+	LWGEOM *l51 = NULL;
 	char *ewkt;
 
 	/* Vertical line with vertices at y integers */
-	l51 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
+	l51 = lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
 
 	/* Clip in the middle, mid-range. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
@@ -562,7 +562,7 @@
 	lwcollection_free(c);
 
 	/* Clip off the top. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 3.5,0 4))");
@@ -570,7 +570,7 @@
 	lwcollection_free(c);
 
 	/* Clip off the bottom. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 2.5))" );
@@ -578,7 +578,7 @@
 	lwcollection_free(c);
 
 	/* Range holds entire object. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))" );
@@ -586,7 +586,7 @@
 	lwcollection_free(c);
 
 	/* Clip on vertices. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1,0 2))" );
@@ -594,7 +594,7 @@
 	lwcollection_free(c);
 
 	/* Clip on vertices off the bottom. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2))" );
@@ -602,7 +602,7 @@
 	lwcollection_free(c);
 
 	/* Clip on top. */
-	c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0);
+	c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(0 0))" );
@@ -610,46 +610,81 @@
 	lwcollection_free(c);
 
 	/* ST_LocateBetweenElevations(ST_GeomFromEWKT('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)'), 1, 2)) */
-	line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
-	c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
+	line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 2.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
 	lwfree(ewkt);
 	lwcollection_free(c);
-	lwline_free(line);
+	lwgeom_free(line);
 
 	/* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 2)) */
-	line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
-	c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
+	line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 2.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("a = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
 	lwfree(ewkt);
 	lwcollection_free(c);
-	lwline_free(line);
+	lwgeom_free(line);
 
 	/* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 1)) */
-	line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
-	c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
+	line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 1.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("b = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
 	lwfree(ewkt);
 	lwcollection_free(c);
-	lwline_free(line);
+	lwgeom_free(line);
 
 	/* ST_LocateBetweenElevations('LINESTRING(1 1 1, 1 2 2)', 1,1) */
-	line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
-	c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
+	line = lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 1.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
 	lwfree(ewkt);
 	lwcollection_free(c);
-	lwline_free(line);
+	lwgeom_free(line);
 
-	lwline_free(l51);
+	lwgeom_free(l51);
+}
 
+static void
+test_lwpoly_clip(void)
+{
+	LWCOLLECTION *c;
+	LWGEOM *g = NULL;
+	char *ewkt;
+
+	g = lwgeom_from_wkt(
+	    "POLYGON ((0.51 -0.25, 1.27 -0.14, 1.27 0.25, 0.6 0.3, 0.7 0.7, 1.2 0.7, 0.8 0.5, 1.3 0.4, 1.2 1.2, 0.5 1.2, 0.5 -0.1, 0.3 -0.1, 0.3 1.3, -0.18 1.25, -0.17 -0.25, 0.51 -0.25))",
+	    LW_PARSER_CHECK_NONE);
+	c = lwgeom_clip_to_ordinate_range(g, 'X', 0.0, 1.0, 0);
+
+	ewkt = lwgeom_to_ewkt((LWGEOM *)c);
+	// printf("c = %s\n", ewkt);
+	CU_ASSERT_STRING_EQUAL(
+	    ewkt,
+	    "MULTIPOLYGON(((0.51 -0.25,1 -0.179078947368,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1.2,0.5 1.2,0.5 -0.1,0.3 -0.1,0.3 1.3,0 1.26875,0 -0.25,0.51 -0.25)))");
+	lwfree(ewkt);
+	lwcollection_free(c);
+	lwgeom_free(g);
+
+	g = lwgeom_from_wkt(
+	    "MULTIPOLYGON(((0.51 -0.25,1 -0.179078947368,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1.2,0.5 1.2,0.5 -0.1,0.3 -0.1,0.3 1.3,0 1.26875,0 -0.25,0.51 -0.25)))",
+	    LW_PARSER_CHECK_NONE);
+	c = lwgeom_clip_to_ordinate_range(g, 'Y', 0.0, 1.0, 0);
+
+	ewkt = lwgeom_to_ewkt((LWGEOM *)c);
+	//printf("c = %s\n", ewkt);
+	CU_ASSERT_STRING_EQUAL(
+	    ewkt,
+	    "MULTIPOLYGON(((1 0,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1,0.5 1,0.5 0,0.3 0,0.3 1,0 1,0 0,1 0)))");
+	lwfree(ewkt);
+	lwcollection_free(c);
+	lwgeom_free(g);
 }
 
 static void test_lwmline_clip(void)
@@ -656,16 +691,16 @@
 {
 	LWCOLLECTION *c;
 	char *ewkt;
-	LWMLINE *mline = NULL;
-	LWLINE *line = NULL;
+	LWGEOM *mline = NULL;
+	LWGEOM *line = NULL;
 
 	/*
 	** Set up the input line. Trivial one-member case.
 	*/
-	mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
+	mline = lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
 
 	/* Clip in the middle, mid-range. */
-	c = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)mline, 'Y', 1.5, 2.5);
+	c = lwgeom_clip_to_ordinate_range(mline, 'Y', 1.5, 2.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
@@ -672,15 +707,15 @@
 	lwfree(ewkt);
 	lwcollection_free(c);
 
-	lwmline_free(mline);
+	lwgeom_free(mline);
 
 	/*
 	** Set up the input line. Two-member case.
 	*/
-	mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((1 0,1 1,1 2,1 3,1 4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
+	mline = lwgeom_from_wkt("MULTILINESTRING((1 0,1 1,1 2,1 3,1 4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
 
 	/* Clip off the top. */
-	c = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)mline, 'Y', 3.5, 5.5);
+	c = lwgeom_clip_to_ordinate_range(mline, 'Y', 3.5, 5.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((1 3.5,1 4),(0 3.5,0 4))");
@@ -687,15 +722,16 @@
 	lwfree(ewkt);
 	lwcollection_free(c);
 
-	lwmline_free(mline);
+	lwgeom_free(mline);
 
 	/*
 	** Set up staggered input line to create multi-type output.
 	*/
-	mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((1 0,1 -1,1 -2,1 -3,1 -4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
+	mline =
+	    lwgeom_from_wkt("MULTILINESTRING((1 0,1 -1,1 -2,1 -3,1 -4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
 
 	/* Clip from 0 upwards.. */
-	c = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)mline, 'Y', 0.0, 2.5);
+	c = lwgeom_clip_to_ordinate_range(mline, 'Y', 0.0, 2.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 0),LINESTRING(0 0,0 1,0 2,0 2.5))");
@@ -702,15 +738,16 @@
 	lwfree(ewkt);
 	lwcollection_free(c);
 
-	lwmline_free(mline);
+	lwgeom_free(mline);
 
 	/*
 	** Set up input line from MAC
 	*/
-	line = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0 0 0,1 1 1 1,2 2 2 2,3 3 3 3,4 4 4 4,3 3 3 5,2 2 2 6,1 1 1 7,0 0 0 8)", LW_PARSER_CHECK_NONE);
+	line = lwgeom_from_wkt("LINESTRING(0 0 0 0,1 1 1 1,2 2 2 2,3 3 3 3,4 4 4 4,3 3 3 5,2 2 2 6,1 1 1 7,0 0 0 8)",
+			       LW_PARSER_CHECK_NONE);
 
 	/* Clip from 3 to 3.5 */
-	c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 3.5);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 3.0, 3.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5))");
@@ -718,7 +755,7 @@
 	lwcollection_free(c);
 
 	/* Clip from 2 to 3.5 */
-	c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.5);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 2.0, 3.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5,2 2 2 6))");
@@ -726,7 +763,7 @@
 	lwcollection_free(c);
 
 	/* Clip from 3 to 4 */
-	c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 4.0);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 3.0, 4.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,4 4 4 4,3 3 3 5))");
@@ -734,7 +771,7 @@
 	lwcollection_free(c);
 
 	/* Clip from 2 to 3 */
-	c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.0);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 2.0, 3.0, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3),(3 3 3 5,2 2 2 6))");
@@ -741,17 +778,13 @@
 	lwfree(ewkt);
 	lwcollection_free(c);
 
-
-	lwline_free(line);
-
+	lwgeom_free(line);
 }
 
-
-
 static void test_lwline_clip_big(void)
 {
 	POINTARRAY *pa = ptarray_construct(1, 0, 3);
-	LWLINE *line = lwline_construct(SRID_UNKNOWN, NULL, pa);
+	LWGEOM *line = (LWGEOM *)lwline_construct(SRID_UNKNOWN, NULL, pa);
 	LWCOLLECTION *c;
 	char *ewkt;
 	POINT4D p;
@@ -771,7 +804,7 @@
 	p.z = 2.0;
 	ptarray_set_point4d(pa, 2, &p);
 
-	c = lwline_clip_to_ordinate_range(line, 'Z', 0.5, 1.5);
+	c = lwgeom_clip_to_ordinate_range(line, 'Z', 0.5, 1.5, 0);
 	ewkt = lwgeom_to_ewkt((LWGEOM*)c);
 	//printf("c = %s\n", ewkt);
 	CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0.5 0.5 0.5,1 1 1,1.5 1.5 1.5))" );
@@ -778,7 +811,7 @@
 
 	lwfree(ewkt);
 	lwcollection_free(c);
-	lwline_free(line);
+	lwgeom_free(line);
 }
 
 static void test_geohash_precision(void)
@@ -1497,6 +1530,7 @@
 	PG_ADD_TEST(suite,test_point_interpolate);
 	PG_ADD_TEST(suite,test_lwline_interpolate_points);
 	PG_ADD_TEST(suite,test_lwline_clip);
+	PG_ADD_TEST(suite, test_lwpoly_clip);
 	PG_ADD_TEST(suite,test_lwline_clip_big);
 	PG_ADD_TEST(suite,test_lwmline_clip);
 	PG_ADD_TEST(suite,test_geohash_point);

Modified: trunk/liblwgeom/liblwgeom_internal.h
===================================================================
--- trunk/liblwgeom/liblwgeom_internal.h	2018-10-25 14:25:19 UTC (rev 16951)
+++ trunk/liblwgeom/liblwgeom_internal.h	2018-10-25 23:16:02 UTC (rev 16952)
@@ -234,26 +234,6 @@
 int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value);
 
 
-/**
-* Clip a line based on the from/to range of one of its ordinates. Use for m- and z- clipping
-*/
-LWCOLLECTION *lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to);
-
-/**
- * Clip collection based on the from/to range of one of its ordinates. Use for m- and z- clipping
- */
-LWCOLLECTION *lwcollection_clip_to_ordinate_range(const LWCOLLECTION *col, char ordinate, double from, double to);
-
-/**
-* Clip a multi-point based on the from/to range of one of its ordinates. Use for m- and z- clipping
-*/
-LWCOLLECTION *lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to);
-
-/**
-* Clip a point based on the from/to range of one of its ordinates. Use for m- and z- clipping
-*/
-LWCOLLECTION *lwpoint_clip_to_ordinate_range(const LWPOINT *mpoint, char ordinate, double from, double to);
-
 /*
 * Geohash
 */

Modified: trunk/liblwgeom/lwlinearreferencing.c
===================================================================
--- trunk/liblwgeom/lwlinearreferencing.c	2018-10-25 14:25:19 UTC (rev 16951)
+++ trunk/liblwgeom/lwlinearreferencing.c	2018-10-25 23:16:02 UTC (rev 16952)
@@ -248,32 +248,28 @@
 * @param ordinate number (1=x, 2=y, 3=z, 4=m)
 * @return d value at that ordinate
 */
-double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
+inline double
+lwpoint_get_ordinate(const POINT4D *p, char ordinate)
 {
-	if ( ! p )
+	if (!p)
 	{
 		lwerror("Null input geometry.");
 		return 0.0;
 	}
 
-	if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
+	switch (ordinate)
 	{
-		lwerror("Cannot extract %c ordinate.", ordinate);
-		return 0.0;
-	}
-
-	if ( ordinate == 'X' )
+	case 'X':
 		return p->x;
-	if ( ordinate == 'Y' )
+	case 'Y':
 		return p->y;
-	if ( ordinate == 'Z' )
+	case 'Z':
 		return p->z;
-	if ( ordinate == 'M' )
+	case 'M':
 		return p->m;
-
-	/* X */
-	return p->x;
-
+	}
+	lwerror("Cannot extract %c ordinate.", ordinate);
+	return 0.0;
 }
 
 /**
@@ -280,22 +276,21 @@
 * Given a point, ordinate number and value, set that ordinate on the
 * point.
 */
-void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
+inline void
+lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
 {
-	if ( ! p )
+	if (!p)
 	{
 		lwerror("Null input geometry.");
 		return;
 	}
 
-	if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
+	if (!(ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M'))
 	{
 		lwerror("Cannot set %c ordinate.", ordinate);
 		return;
 	}
 
-	LWDEBUGF(4, "    setting ordinate %c to %g", ordinate, value);
-
 	switch ( ordinate )
 	{
 	case 'X':
@@ -318,7 +313,14 @@
 * generate a new point that is proportionally between the input points,
 * using the values in the provided dimension as the scaling factors.
 */
-int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
+inline int
+point_interpolate(const POINT4D *p1,
+		  const POINT4D *p2,
+		  POINT4D *p,
+		  int hasz,
+		  int hasm,
+		  char ordinate,
+		  double interpolation_value)
 {
 	static char* dims = "XYZM";
 	double p1_value = lwpoint_get_ordinate(p1, ordinate);
@@ -326,34 +328,44 @@
 	double proportion;
 	int i = 0;
 
-	if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
+#if PARANOIA_LEVEL > 0
+	if (!(ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M'))
 	{
-		lwerror("Cannot set %c ordinate.", ordinate);
-		return 0;
+		lwerror("Cannot interpolate over %c ordinate.", ordinate);
+		return LW_FAILURE;
 	}
 
-	if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
-	        FP_MAX(p1_value, p2_value) < interpolation_value )
+	if (FP_MIN(p1_value, p2_value) > interpolation_value || FP_MAX(p1_value, p2_value) < interpolation_value)
 	{
-		lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
-		return 0;
+		lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).",
+			interpolation_value,
+			p1_value,
+			p2_value);
+		return LW_FAILURE;
 	}
+#endif
 
 	proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
 
-	for ( i = 0; i < 4; i++ )
+	for (i = 0; i < 4; i++)
 	{
 		double newordinate = 0.0;
-		if ( dims[i] == 'Z' && ! hasz ) continue;
-		if ( dims[i] == 'M' && ! hasm ) continue;
-		p1_value = lwpoint_get_ordinate(p1, dims[i]);
-		p2_value = lwpoint_get_ordinate(p2, dims[i]);
-		newordinate = p1_value + proportion * (p2_value - p1_value);
-		lwpoint_set_ordinate(p, dims[i], newordinate);
-		LWDEBUGF(4, "   clip ordinate(%c) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", dims[i], p1_value, p2_value, proportion, newordinate );
+		if (dims[i] == 'Z' && !hasz)
+			continue;
+		if (dims[i] == 'M' && !hasm)
+			continue;
+		if (dims[i] == ordinate)
+			lwpoint_set_ordinate(p, dims[i], interpolation_value);
+		else
+		{
+			p1_value = lwpoint_get_ordinate(p1, dims[i]);
+			p2_value = lwpoint_get_ordinate(p2, dims[i]);
+			newordinate = p1_value + proportion * (p2_value - p1_value);
+			lwpoint_set_ordinate(p, dims[i], newordinate);
+		}
 	}
 
-	return 1;
+	return LW_SUCCESS;
 }
 
 
@@ -360,7 +372,7 @@
 /**
 * Clip an input POINT between two values, on any ordinate input.
 */
-LWCOLLECTION*
+static inline LWCOLLECTION *
 lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
 {
 	LWCOLLECTION *lwgeom_out = NULL;
@@ -368,18 +380,6 @@
 	POINT4D p4d;
 	double ordinate_value;
 
-	/* Nothing to do with NULL */
-	if ( ! point )
-		lwerror("Null input geometry.");
-
-	/* Ensure 'from' is less than 'to'. */
-	if ( to < from )
-	{
-		double t = from;
-		from = to;
-		to = t;
-	}
-
 	/* Read Z/M info */
 	hasz = lwgeom_has_z(lwpoint_as_lwgeom(point));
 	hasm = lwgeom_has_m(lwpoint_as_lwgeom(point));
@@ -396,21 +396,13 @@
 		lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
 	}
 
-	/* Set the bbox, if necessary */
-	if ( lwgeom_out->bbox )
-	{
-		lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
-	}
-
 	return lwgeom_out;
 }
 
-
-
 /**
 * Clip an input MULTIPOINT between two values, on any ordinate input.
 */
-LWCOLLECTION*
+static inline LWCOLLECTION *
 lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
 {
 	LWCOLLECTION *lwgeom_out = NULL;
@@ -417,18 +409,6 @@
 	char hasz, hasm;
 	uint32_t i;
 
-	/* Nothing to do with NULL */
-	if ( ! mpoint )
-		lwerror("Null input geometry.");
-
-	/* Ensure 'from' is less than 'to'. */
-	if ( to < from )
-	{
-		double t = from;
-		from = to;
-		to = t;
-	}
-
 	/* Read Z/M info */
 	hasz = lwgeom_has_z(lwmpoint_as_lwgeom(mpoint));
 	hasm = lwgeom_has_m(lwmpoint_as_lwgeom(mpoint));
@@ -453,62 +433,112 @@
 	}
 
 	/* Set the bbox, if necessary */
-	if ( lwgeom_out->bbox )
-	{
+	if (mpoint->bbox)
 		lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
-	}
 
 	return lwgeom_out;
 }
 
-/**
- * Clip an input COLLECTION between two values, on any ordinate input.
- */
-LWCOLLECTION *
-lwcollection_clip_to_ordinate_range(const LWCOLLECTION *icol, char ordinate, double from, double to)
+static inline POINTARRAY *
+ptarray_clamp_to_ordinate_range(const POINTARRAY *ipa, char ordinate, double from, double to)
 {
-	LWCOLLECTION *lwgeom_out = NULL;
+	POINT4D p1, p2;
+	POINTARRAY *opa;
+	double ovp1, ovp2;
+	POINT4D *t;
+	int8_t p1out, p2out; /* -1 - smaller than from, 0 - in range, 1 - larger than to */
+	uint32_t i;
+	uint8_t hasz = FLAGS_GET_Z(ipa->flags);
+	uint8_t hasm = FLAGS_GET_M(ipa->flags);
 
-	if (!icol)
+	assert(from <= to);
+
+	t = lwalloc(sizeof(POINT4D));
+
+	/* Initial storage */
+	opa = ptarray_construct_empty(hasz, hasm, ipa->npoints);
+
+	/* Add first point */
+	getPoint4d_p(ipa, 0, &p1);
+	ovp1 = lwpoint_get_ordinate(&p1, ordinate);
+
+	p1out = (ovp1 < from) ? -1 : ((ovp1 > to) ? 1 : 0);
+
+	if (from <= ovp1 && ovp1 <= to)
+		ptarray_append_point(opa, &p1, LW_FALSE);
+
+	/* Loop on all other input points */
+	for (i = 1; i < ipa->npoints; i++)
 	{
-		lwerror("Null input geometry.");
-		return NULL;
-	}
+		getPoint4d_p(ipa, i, &p2);
+		ovp2 = lwpoint_get_ordinate(&p2, ordinate);
+		p2out = (ovp2 < from) ? -1 : ((ovp2 > to) ? 1 : 0);
 
-	if (icol->ngeoms == 1)
-		lwgeom_out = lwgeom_clip_to_ordinate_range(icol->geoms[0], ordinate, from, to, 0);
-	else
-	{
-		LWCOLLECTION *col;
-		char hasz = lwgeom_has_z(lwcollection_as_lwgeom(icol));
-		char hasm = lwgeom_has_m(lwcollection_as_lwgeom(icol));
-		uint32_t i;
-		lwgeom_out = lwcollection_construct_empty(icol->type, icol->srid, hasz, hasm);
-		FLAGS_SET_Z(lwgeom_out->flags, hasz);
-		FLAGS_SET_M(lwgeom_out->flags, hasm);
-		for (i = 0; i < icol->ngeoms; i++)
+		if (p1out == 0 && p2out == 0) /* both visible */
 		{
-			col = lwgeom_clip_to_ordinate_range(icol->geoms[i], ordinate, from, to, 0);
-			if (col)
-			{
-				if (col->type != icol->type)
-					lwgeom_out->type = COLLECTIONTYPE;
-				lwgeom_out = lwcollection_concat_in_place(lwgeom_out, col);
-				lwcollection_release(col);
-			}
+			ptarray_append_point(opa, &p2, LW_FALSE);
 		}
-		if (lwgeom_out->bbox)
-			lwgeom_refresh_bbox((LWGEOM *)lwgeom_out);
+		else if (p1out == p2out && p1out != 0) /* both invisible */
+		{
+			/* skip */
+		}
+		else if (p1out == -1 && p2out == 0)
+		{
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
+			ptarray_append_point(opa, t, LW_FALSE);
+			ptarray_append_point(opa, &p2, LW_FALSE);
+		}
+		else if (p1out == -1 && p2out == 1)
+		{
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
+			ptarray_append_point(opa, t, LW_FALSE);
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
+			ptarray_append_point(opa, t, LW_FALSE);
+		}
+		else if (p1out == 0 && p2out == -1)
+		{
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
+			ptarray_append_point(opa, t, LW_FALSE);
+		}
+		else if (p1out == 0 && p2out == 1)
+		{
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
+			ptarray_append_point(opa, t, LW_FALSE);
+		}
+		else if (p1out == 1 && p2out == -1)
+		{
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
+			ptarray_append_point(opa, t, LW_FALSE);
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
+			ptarray_append_point(opa, t, LW_FALSE);
+		}
+		else if (p1out == 1 && p2out == 0)
+		{
+			point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
+			ptarray_append_point(opa, t, LW_FALSE);
+			ptarray_append_point(opa, &p2, LW_FALSE);
+		}
+
+		p1 = p2;
+		p1out = p2out;
+		LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
 	}
-	return lwgeom_out;
+
+	if (ptarray_is_closed(ipa) && opa->npoints > 0)
+	{
+		getPoint4d_p(opa, 0, &p1);
+		ptarray_append_point(opa, &p1, LW_FALSE);
+	}
+	lwfree(t);
+
+	return opa;
 }
 
-
 /**
 * Take in a LINESTRING and return a MULTILINESTRING of those portions of the
 * LINESTRING between the from/to range for the specified ordinate (XYZM)
 */
-LWCOLLECTION*
+static inline LWCOLLECTION *
 lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
 {
 
@@ -535,21 +565,6 @@
 	hasm = lwgeom_has_m(lwline_as_lwgeom(line));
 	dims = FLAGS_NDIMS(line->flags);
 
-	/* Ensure 'from' is less than 'to'. */
-	if ( to < from )
-	{
-		double t = from;
-		from = to;
-		to = t;
-	}
-
-#if POSTGIS_DEBUG_LEVEL >= 4
-	LWDEBUGF(4, "from = %g, to = %g, ordinate = %c", from, to, ordinate);
-	geom_ewkt = lwgeom_to_ewkt((LWGEOM*)line);
-	LWDEBUGF(4, "%s", geom_ewkt);
-	lwfree(geom_ewkt);
-#endif
-
 	/* Asking for an ordinate we don't have. Error. */
 	if ( (ordinate == 'Z' && ! hasz) || (ordinate == 'M' && ! hasm) )
 	{
@@ -730,13 +745,94 @@
 	lwfree(q);
 	lwfree(r);
 
-	if ( lwgeom_out->bbox && lwgeom_out->ngeoms > 0 )
+	if (line->bbox && lwgeom_out->ngeoms > 0)
+		lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
+
+	return lwgeom_out;
+}
+
+/**
+ * Clip an input LWPOLY between two values, on any ordinate input.
+ */
+static inline LWCOLLECTION *
+lwpoly_clip_to_ordinate_range(const LWPOLY *poly, char ordinate, double from, double to)
+{
+	LWCOLLECTION *lwgeom_out = NULL;
+	uint32_t i, nrings;
+	char hasz = FLAGS_GET_Z(poly->flags), hasm = FLAGS_GET_M(poly->flags);
+	LWPOLY *poly_res = lwpoly_construct_empty(poly->srid, hasz, hasm);
+
+	assert(poly);
+	lwgeom_out = lwcollection_construct_empty(MULTIPOLYGONTYPE, poly->srid, hasz, hasm);
+
+	nrings = poly->nrings;
+	for (i = 0; i < nrings; i++)
 	{
-		lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
+		/* Ret number of points */
+		POINTARRAY *pa = ptarray_clamp_to_ordinate_range(poly->rings[i], ordinate, from, to);
+
+		if (pa->npoints >= 4)
+			lwpoly_add_ring(poly_res, pa);
+		else
+		{
+			ptarray_free(pa);
+			if (i == 0)
+				break;
+		}
 	}
+	lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, (LWGEOM *)poly_res);
 
 	return lwgeom_out;
+}
 
+/**
+ * Clip an input COLLECTION between two values, on any ordinate input.
+ */
+static inline LWCOLLECTION *
+lwcollection_clip_to_ordinate_range(const LWCOLLECTION *icol, char ordinate, double from, double to)
+{
+	LWCOLLECTION *lwgeom_out = NULL;
+
+	if (!icol)
+	{
+		lwerror("Null input geometry.");
+		return NULL;
+	}
+
+	/* Ensure 'from' is less than 'to'. */
+	if (to < from)
+	{
+		double t = from;
+		from = to;
+		to = t;
+	}
+
+	if (icol->ngeoms == 1)
+		lwgeom_out = lwgeom_clip_to_ordinate_range(icol->geoms[0], ordinate, from, to, 0);
+	else
+	{
+		LWCOLLECTION *col;
+		char hasz = lwgeom_has_z(lwcollection_as_lwgeom(icol));
+		char hasm = lwgeom_has_m(lwcollection_as_lwgeom(icol));
+		uint32_t i;
+		lwgeom_out = lwcollection_construct_empty(icol->type, icol->srid, hasz, hasm);
+		FLAGS_SET_Z(lwgeom_out->flags, hasz);
+		FLAGS_SET_M(lwgeom_out->flags, hasm);
+		for (i = 0; i < icol->ngeoms; i++)
+		{
+			col = lwgeom_clip_to_ordinate_range(icol->geoms[i], ordinate, from, to, 0);
+			if (col)
+			{
+				if (col->type != icol->type)
+					lwgeom_out->type = COLLECTIONTYPE;
+				lwgeom_out = lwcollection_concat_in_place(lwgeom_out, col);
+				lwcollection_release(col);
+			}
+		}
+		if (icol->bbox)
+			lwgeom_refresh_bbox((LWGEOM *)lwgeom_out);
+	}
+	return lwgeom_out;
 }
 
 LWCOLLECTION*
@@ -760,6 +856,9 @@
 	case POINTTYPE:
 		out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to);
 		break;
+	case POLYGONTYPE:
+		out_col = lwpoly_clip_to_ordinate_range((LWPOLY *)lwin, ordinate, from, to);
+		break;
 	// case TRIANGLETYPE:
 	// 	out_col = lwtriangle_clip_to_ordinate_range((LWTRIANGLE*)lwin, ordinate, from, to);
 	// 	break;
@@ -775,7 +874,7 @@
 	}
 
 	/* Stop if result is NULL */
-	if ( out_col == NULL )
+	if (!out_col)
 		lwerror("lwgeom_clip_to_ordinate_range clipping routine returned NULL");
 
 	/* Return if we aren't going to offset the result */



More information about the postgis-tickets mailing list