Commit b79db62a authored by Darafei Praliaskouski's avatar Darafei Praliaskouski

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
parent 8b396228
Pipeline #28510126 passed with stage
in 23 minutes and 37 seconds
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
......
......@@ -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
......
......@@ -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);
......
......@@ -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);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment