[postgis-tickets] r17097 - ST_3DLineInterpolatePoint.
Darafei
komzpa at gmail.com
Mon Dec 3 12:16:41 PST 2018
Author: komzpa
Date: 2018-12-03 12:16:40 -0800 (Mon, 03 Dec 2018)
New Revision: 17097
Modified:
trunk/NEWS
trunk/doc/reference_lrs.xml
trunk/liblwgeom/cunit/cu_algorithm.c
trunk/liblwgeom/liblwgeom.h.in
trunk/liblwgeom/lwline.c
trunk/postgis/lwgeom_functions_analytic.c
trunk/postgis/postgis.sql.in
trunk/regress/core/regress_lrs.sql
trunk/regress/core/regress_lrs_expected
Log:
ST_3DLineInterpolatePoint.
Patch by Julien Cabieces and Vincent Mora.
Closes #2645
Closes #4171
Closes https://github.com/postgis/postgis/pull/344
Closes https://github.com/postgis/postgis/pull/346
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/NEWS 2018-12-03 20:16:40 UTC (rev 17097)
@@ -17,6 +17,7 @@
- #4230, SP-GiST and GiST support for ND box operators overlaps, contains,
within, equals (Esteban Zimányi and Arthur Lesuisse from Université
Libre de Bruxelles (ULB), Darafei Praliaskouski)
+ - #4171, ST_3DLineInterpolatePoint (Julien Cabieces, Vincent Mora)
* Enhancements and fixes *
- #4153, ST_Segmentize now splits segments proportionally (Darafei
Modified: trunk/doc/reference_lrs.xml
===================================================================
--- trunk/doc/reference_lrs.xml 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/doc/reference_lrs.xml 2018-12-03 20:16:40 UTC (rev 17097)
@@ -88,11 +88,75 @@
<xref linkend="ST_AsEWKT" />,
<xref linkend="ST_Length" />,
<xref linkend="ST_LineInterpolatePoints" />
+ <xref linkend="ST_3DLineInterpolatePoint" />
<xref linkend="ST_LineLocatePoint" />
O</para>
</refsection>
</refentry>
+ <refentry id="ST_3DLineInterpolatePoint">
+ <refnamediv>
+ <refname>ST_3DLineInterpolatePoint</refname>
+
+ <refpurpose>Returns a point interpolated along a line in 3D. Second argument is a float8 between 0 and 1
+ representing fraction of total length of linestring the point has to be located.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_LineInterpolatePoint</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>a_linestring</parameter></paramdef>
+ <paramdef><type>float8 </type> <parameter>a_fraction</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Returns a point interpolated along a line. First argument
+ must be a LINESTRING. Second argument is a float8 between 0 and 1
+ representing fraction of total linestring length the point has to be located.</para>
+
+ <note>
+ <para><xref linkend="ST_LineInterpolatePoint" /> computes resulting point in 2D and then interpolates
+ value for Z and M, while ST_3DLineInterpolatePoint computes directly point in 3D and only M value
+ is interpolated then.</para>
+ </note>
+
+ <para>Availability: 3.0.0</para>
+ </refsection>
+
+
+ <refsection>
+ <title>Examples</title>
+
+ <para>Return point 20% along 3D line</para>
+ <programlisting>
+SELECT ST_AsEWKT(ST_3DLineInterpolatePoint(the_line, 0.20))
+ FROM (SELECT ST_GeomFromEWKT('LINESTRING(25 50 70, 100 125 90, 150 190 200)') as the_line) As foo;
+ st_asewkt
+----------------
+ POINT(59.0675892910822 84.0675892910822 79.0846904776219)
+</programlisting>
+ </refsection>
+
+ <!-- Optionally add a "See Also" section -->
+ <refsection>
+ <title>See Also</title>
+
+ <para>
+ <xref linkend="ST_AsText" />,
+ <xref linkend="ST_AsEWKT" />,
+ <xref linkend="ST_Length" />,
+ <xref linkend="ST_LineInterpolatePoint" />
+ <xref linkend="ST_LineInterpolatePoints" />
+ <xref linkend="ST_LineLocatePoint" />
+ </para>
+ </refsection>
+ </refentry>
+
<refentry id="ST_LineInterpolatePoints">
<refnamediv>
<refname>ST_LineInterpolatePoints</refname>
Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/liblwgeom/cunit/cu_algorithm.c 2018-12-03 20:16:40 UTC (rev 17097)
@@ -543,6 +543,59 @@
lwgeom_free(lwline_as_lwgeom(empty_line));
}
+static void
+test_lwline_interpolate_point_3d(void)
+{
+ LWLINE *line;
+ POINT4D point;
+ LWPOINT *pt;
+
+ /* Empty line -> Empty point*/
+ line = lwline_construct_empty(4326, LW_TRUE, LW_FALSE);
+ pt = lwline_interpolate_point_3d(line, 0.5);
+ CU_ASSERT(lwpoint_is_empty(pt));
+ CU_ASSERT(lwgeom_has_z(lwpoint_as_lwgeom(pt)));
+ CU_ASSERT_FALSE(lwgeom_has_m(lwpoint_as_lwgeom(pt)));
+ lwpoint_free(pt);
+ lwline_free(line);
+
+ line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING Z (0 0 0, 1 1 1, 2 2 2)", LW_PARSER_CHECK_NONE));
+
+ /* distance = 0 -> first point */
+ pt = lwline_interpolate_point_3d(line, 0);
+ lwpoint_getPoint4d_p(pt, &point);
+ CU_ASSERT_DOUBLE_EQUAL(point.x, 0, 0.0001);
+ CU_ASSERT_DOUBLE_EQUAL(point.y, 0, 0.001);
+ CU_ASSERT_DOUBLE_EQUAL(point.z, 0, 0.001);
+ lwpoint_free(pt);
+
+ /* distance = 1 -> last point */
+ pt = lwline_interpolate_point_3d(line, 1);
+ lwpoint_getPoint4d_p(pt, &point);
+ CU_ASSERT_DOUBLE_EQUAL(point.x, 2, 0.0001);
+ CU_ASSERT_DOUBLE_EQUAL(point.y, 2, 0.001);
+ CU_ASSERT_DOUBLE_EQUAL(point.z, 2, 0.001);
+ lwpoint_free(pt);
+
+ /* simple where distance 50% -> second point */
+ pt = lwline_interpolate_point_3d(line, 0.5);
+ lwpoint_getPoint4d_p(pt, &point);
+ CU_ASSERT_DOUBLE_EQUAL(point.x, 1, 0.0001);
+ CU_ASSERT_DOUBLE_EQUAL(point.y, 1, 0.001);
+ CU_ASSERT_DOUBLE_EQUAL(point.z, 1, 0.001);
+ lwpoint_free(pt);
+
+ /* simple where distance 80% -> between second and last point */
+ pt = lwline_interpolate_point_3d(line, 0.8);
+ lwpoint_getPoint4d_p(pt, &point);
+ CU_ASSERT_DOUBLE_EQUAL(point.x, 1.6, 0.0001);
+ CU_ASSERT_DOUBLE_EQUAL(point.y, 1.6, 0.001);
+ CU_ASSERT_DOUBLE_EQUAL(point.z, 1.6, 0.001);
+ lwpoint_free(pt);
+
+ lwline_free(line);
+}
+
static void test_lwline_clip(void)
{
LWCOLLECTION *c;
@@ -1564,6 +1617,7 @@
PG_ADD_TEST(suite,test_lwpoint_get_ordinate);
PG_ADD_TEST(suite,test_point_interpolate);
PG_ADD_TEST(suite,test_lwline_interpolate_points);
+ PG_ADD_TEST(suite, test_lwline_interpolate_point_3d);
PG_ADD_TEST(suite,test_lwline_clip);
PG_ADD_TEST(suite, test_lwpoly_clip);
PG_ADD_TEST(suite, test_lwtriangle_clip);
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/liblwgeom/liblwgeom.h.in 2018-12-03 20:16:40 UTC (rev 17097)
@@ -1002,6 +1002,11 @@
*/
extern POINTARRAY* lwline_interpolate_points(const LWLINE *line, double length_fraction, char repeat);
+/**
+ * Interpolate one point along a line in 3D
+ */
+extern LWPOINT* lwline_interpolate_point_3d(const LWLINE *line, double distance);
+
/******************************************************************
* LWPOLY functions
******************************************************************/
Modified: trunk/liblwgeom/lwline.c
===================================================================
--- trunk/liblwgeom/lwline.c 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/liblwgeom/lwline.c 2018-12-03 20:16:40 UTC (rev 17097)
@@ -607,3 +607,69 @@
return opa;
}
+extern LWPOINT *
+lwline_interpolate_point_3d(const LWLINE *line, double distance)
+{
+ double length, slength, tlength;
+ POINTARRAY *ipa;
+ POINT4D pt;
+ int nsegs, i;
+ LWGEOM *geom = lwline_as_lwgeom(line);
+ int has_z = lwgeom_has_z(geom);
+ int has_m = lwgeom_has_m(geom);
+ ipa = line->points;
+
+ /* Empty.InterpolatePoint == Point Empty */
+ if (lwline_is_empty(line))
+ {
+ return lwpoint_construct_empty(line->srid, has_z, has_m);
+ }
+
+ /* If distance is one of the two extremes, return the point on that
+ * end rather than doing any expensive computations
+ */
+ if (distance == 0.0 || distance == 1.0)
+ {
+ if (distance == 0.0)
+ getPoint4d_p(ipa, 0, &pt);
+ else
+ getPoint4d_p(ipa, ipa->npoints - 1, &pt);
+
+ return lwpoint_make(line->srid, has_z, has_m, &pt);
+ }
+
+ /* Interpolate a point on the line */
+ nsegs = ipa->npoints - 1;
+ length = ptarray_length(ipa);
+ tlength = 0;
+ for (i = 0; i < nsegs; i++)
+ {
+ POINT4D p1, p2;
+ POINT4D *p1ptr = &p1, *p2ptr = &p2; /* don't break
+ * strict-aliasing rules
+ */
+
+ getPoint4d_p(ipa, i, &p1);
+ getPoint4d_p(ipa, i + 1, &p2);
+
+ /* Find the relative length of this segment */
+ slength = distance3d_pt_pt((POINT3D *)p1ptr, (POINT3D *)p2ptr) / length;
+
+ /* If our target distance is before the total length we've seen
+ * so far. create a new point some distance down the current
+ * segment.
+ */
+ if (distance < tlength + slength)
+ {
+ double dseg = (distance - tlength) / slength;
+ interpolate_point4d(&p1, &p2, &pt, dseg);
+ return lwpoint_make(line->srid, has_z, has_m, &pt);
+ }
+ tlength += slength;
+ }
+
+ /* Return the last point on the line. This shouldn't happen, but
+ * could if there's some floating point rounding errors. */
+ getPoint4d_p(ipa, ipa->npoints - 1, &pt);
+ return lwpoint_make(line->srid, has_z, has_m, &pt);
+}
Modified: trunk/postgis/lwgeom_functions_analytic.c
===================================================================
--- trunk/postgis/lwgeom_functions_analytic.c 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/postgis/lwgeom_functions_analytic.c 2018-12-03 20:16:40 UTC (rev 17097)
@@ -223,7 +223,51 @@
PG_RETURN_POINTER(result);
}
+
/***********************************************************************
+ * Interpolate a point along a line 3D version
+ * --vincent.mora at oslandia.com;
+ ***********************************************************************/
+
+Datum ST_3DLineInterpolatePoint(PG_FUNCTION_ARGS);
+
+PG_FUNCTION_INFO_V1(ST_3DLineInterpolatePoint);
+Datum ST_3DLineInterpolatePoint(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gser = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *result;
+ double distance = PG_GETARG_FLOAT8(1);
+ LWLINE *line;
+ LWGEOM *geom;
+ LWPOINT *point;
+
+ if (distance < 0 || distance > 1)
+ {
+ elog(ERROR, "line_interpolate_point: 2nd arg isn't within [0,1]");
+ PG_RETURN_NULL();
+ }
+
+ if (gserialized_get_type(gser) != LINETYPE)
+ {
+ elog(ERROR, "line_interpolate_point: 1st arg isn't a line");
+ PG_RETURN_NULL();
+ }
+
+ geom = lwgeom_from_gserialized(gser);
+ line = lwgeom_as_lwline(geom);
+
+ point = lwline_interpolate_point_3d(line, distance);
+
+ lwgeom_free(geom);
+ PG_FREE_IF_COPY(gser, 0);
+
+ result = geometry_serialize(lwpoint_as_lwgeom(point));
+ lwpoint_free(point);
+
+ PG_RETURN_POINTER(result);
+}
+
+/***********************************************************************
* --jsunday at rochgrp.com;
***********************************************************************/
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/postgis/postgis.sql.in 2018-12-03 20:16:40 UTC (rev 17097)
@@ -6243,7 +6243,13 @@
GRANT SELECT ON TABLE geometry_columns TO public;
GRANT SELECT ON TABLE spatial_ref_sys TO public;
+-- Availability: 3.0.0
+CREATE OR REPLACE FUNCTION ST_3DLineInterpolatePoint(geometry, float8)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'ST_3DLineInterpolatePoint'
+ LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL;
+
-- moved to separate file cause its invovled
#include "postgis_spgist.sql.in"
Modified: trunk/regress/core/regress_lrs.sql
===================================================================
--- trunk/regress/core/regress_lrs.sql 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/regress/core/regress_lrs.sql 2018-12-03 20:16:40 UTC (rev 17097)
@@ -73,6 +73,13 @@
select 'line_interpolate_point', ST_AsText(ST_LineInterpolatePoint('LINESTRING(0 0 10, 1 1 5)', 0.5));
--
+--- ST_3DLineInterpolatePoint
+--
+
+select 'line_interpolate_point_3d', ST_AsText(ST_3DLineInterpolatePoint('LINESTRING(0 0 0, 0 0 1)', 0));
+select 'line_interpolate_point_3d', ST_AsText(ST_3DLineInterpolatePoint('LINESTRING(0 0 0, 0 0 1)', 0.5));
+
+--
--- ST_LineInterpolatePoints
--
Modified: trunk/regress/core/regress_lrs_expected
===================================================================
--- trunk/regress/core/regress_lrs_expected 2018-12-03 20:08:58 UTC (rev 17096)
+++ trunk/regress/core/regress_lrs_expected 2018-12-03 20:16:40 UTC (rev 17097)
@@ -39,6 +39,8 @@
line_substring_12|POINT Z (0.5 0.5 7.5)
line_interpolate_point|POINT(0 0)
line_interpolate_point|POINT Z (0.5 0.5 7.5)
+line_interpolate_point_3d|POINT Z (0 0 0)
+line_interpolate_point_3d|POINT Z (0 0 0.5)
line_interpolate_points|POINT(0.7 0.7)
line_interpolate_points|MULTIPOINT(0.3 0.3,0.6 0.6,0.9 0.9)
addMeasure1|LINESTRING M (0 0 10,2 0 15,4 0 20)
More information about the postgis-tickets
mailing list