ST_Segmentize: split geometry proportionally.

This helps avoid situations when 52 meter long segment is split into 50 and 2 meter long ones,
and split into 26 and 26 meter long ones. It is crucial for good interpolating M and Z values.

Thanks Andrew Shadoura for supporting me on Patreon: https://www.patreon.com/komzpa

Closes https://github.com/postgis/postgis/pull/287
Closes #4153

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

diff --git a/NEWS b/NEWS
index ab41cc2c6c471a36f5d56411d2148d0904fd625e..4465dfdc46da211ca7c2c6647acf5ffa2ba725c1 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,10 @@
+PostGIS 3.0.0
+2019/xx/xx
+* Enhancements *
+  - #4153, ST_Segmentize now splits segments proportionally (Darafei
+    Praliaskouski).
+
+
 PostGIS 2.5.0rc1
 2018/08/19
 New since PostGIS 2.5.0beta2

diff --git a/doc/reference_editor.xml b/doc/reference_editor.xml
index 41606ca302447f9a192fd41d9f1bf3e652889ccf..f7d48e73752e5319cf360a781c031e03ba80d0de 100644 (file)
--- a/doc/reference_editor.xml
+++ b/doc/reference_editor.xml
@@ -1556,6 +1556,7 @@ SELECT ST_AsText(ST_Scale('LINESTRING(1 1, 2 2)', 'POINT(2 2)', 'POINT(1 1)'::ge
                        given <varname>max_segment_length</varname>.  Distance computation is performed in 2d
                        only. For geometry, length units are in units of spatial reference.  For geography, units are in meters.</para>
                        <para>Availability: 1.2.2</para>
+                       <para>Enhanced: 3.0.0 Segmentize geometry now uses equal length segments</para>
                        <para>Enhanced: 2.3.0 Segmentize geography now uses equal length segments</para>
                        <para>Enhanced: 2.1.0 support for geography was introduced.</para>
                        <para>Changed: 2.1.0 As a result of the introduction of geography support: The construct <code>SELECT ST_Segmentize('LINESTRING(1 2, 3 4)',0.5);</code> will result in ambiguous function error.  You need to have properly typed object e.g. a geometry/geography column, use ST_GeomFromText, ST_GeogFromText or

diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c
index 6ca9d0146b3f3f79a184a2ea9b24352e2909a62c..569a7e9c4d002b18c87a00584bfd58a1447f4a3e 100644 (file)
--- a/liblwgeom/cunit/cu_measures.c
+++ b/liblwgeom/cunit/cu_measures.c
@@ -518,6 +518,15 @@ test_lwgeom_segmentize2d(void)
        lwgeom_free(linein);
        lwgeom_free(lineout);
 
+       /* test that segmentize is proportional - request every 6, get every 5 */
+       linein = lwgeom_from_wkt("LINESTRING(0 0, 20 0)", LW_PARSER_CHECK_NONE);
+       lineout = lwgeom_segmentize2d(linein, 6);
+       strout = lwgeom_to_ewkt(lineout);
+       ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0,15 0,20 0)");
+       lwfree(strout);
+       lwgeom_free(linein);
+       lwgeom_free(lineout);
+
        /* test interruption */
 
        linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);

diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index a2964426306024a7e576d6e3efef1846ba48f524..1e51119cd9e4b1357da933bba67148db85360234 100644 (file)
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -402,22 +402,21 @@ ptarray_swap_ordinates(POINTARRAY *pa, LWORD o1, LWORD o2)
        }
 }
 
-
 /**
  * @brief Returns a modified #POINTARRAY so that no segment is
  *             longer than the given distance (computed using 2d).
  *
  * Every input point is kept.
- * Z and M values for added points (if needed) are set to 0.
+ * Z and M values for added points (if needed) are set proportionally.
  */
 POINTARRAY *
 ptarray_segmentize2d(const POINTARRAY *ipa, double dist)
 {
-       double  segdist;
+       double segdist;
        POINT4D p1, p2;
        POINT4D pbuf;
        POINTARRAY *opa;
-       uint32_t ipoff=0; /* input point offset */
+       uint32_t i, j, nseg;
        int hasz = FLAGS_GET_Z(ipa->flags);
        int hasm = FLAGS_GET_M(ipa->flags);
 
@@ -427,12 +426,11 @@ ptarray_segmentize2d(const POINTARRAY *ipa, double dist)
        opa = ptarray_construct_empty(hasz, hasm, ipa->npoints);
 
        /* Add first point */
-       getPoint4d_p(ipa, ipoff, &p1);
+       getPoint4d_p(ipa, 0, &p1);
        ptarray_append_point(opa, &p1, LW_FALSE);
 
-       ipoff++;
-
-       while (ipoff<ipa->npoints)
+       /* Loop on all other input points */
+       for (i = 1; i < ipa->npoints; i++)
        {
                /*
                 * We use these pointers to avoid
@@ -442,32 +440,29 @@ ptarray_segmentize2d(const POINTARRAY *ipa, double dist)
                 * It looks that casting a variable address (also
                 * referred to as "type-punned pointer")
                 * breaks those "strict" rules.
-                *
                 */
                POINT4D *p1ptr=&p1, *p2ptr=&p2;
 
-               getPoint4d_p(ipa, ipoff, &p2);
+               getPoint4d_p(ipa, i, &p2);
 
                segdist = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
+               /* Split input segment into shorter even chunks */
+               nseg = ceil(segdist / dist);
 
-               if (segdist > dist) /* add an intermediate point */
+               for (j = 1; j < nseg; j++)
                {
-                       pbuf.x = p1.x + (p2.x-p1.x)/segdist * dist;
-                       pbuf.y = p1.y + (p2.y-p1.y)/segdist * dist;
-                       if( hasz )
-                               pbuf.z = p1.z + (p2.z-p1.z)/segdist * dist;
-                       if( hasm )
-                               pbuf.m = p1.m + (p2.m-p1.m)/segdist * dist;
+                       pbuf.x = p1.x + (p2.x - p1.x) * j / nseg;
+                       pbuf.y = p1.y + (p2.y - p1.y) * j / nseg;
+                       if (hasz)
+                               pbuf.z = p1.z + (p2.z - p1.z) * j / nseg;
+                       if (hasm)
+                               pbuf.m = p1.m + (p2.m - p1.m) * j / nseg;
                        ptarray_append_point(opa, &pbuf, LW_FALSE);
-                       p1 = pbuf;
-               }
-               else /* copy second point */
-               {
-                       ptarray_append_point(opa, &p2, (ipa->npoints==2)?LW_TRUE:LW_FALSE);
-                       p1 = p2;
-                       ipoff++;
+                       LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
                }
 
+               ptarray_append_point(opa, &p2, (ipa->npoints == 2) ? LW_TRUE : LW_FALSE);
+               p1 = p2;
                LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
        }