postgis · Nov 30, 2018 · Nov 30, 2018 · Dec 3, 2018 · Dec 3, 2018 · Dec 3, 2018
diff --git a/doc/reference_lrs.xml b/doc/reference_lrs.xml
@@ -88,11 +88,75 @@ FROM (SELECT ST_GeomFromText('LINESTRING(1 2, 4 5, 6 7)') As the_line) As foo;
 				<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>

diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
@@ -543,6 +543,59 @@ static void test_lwline_interpolate_points(void)
 	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 @@ void algorithms_suite_setup(void)
 	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);

diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
@@ -1002,6 +1002,11 @@ extern int lwline_add_lwpoint(LWLINE *line, LWPOINT *point, uint32_t where);
  */
 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
  ******************************************************************/

diff --git a/liblwgeom/lwline.c b/liblwgeom/lwline.c
@@ -607,3 +607,69 @@ POINTARRAY* lwline_interpolate_points(const LWLINE *line, double length_fraction
     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);
+}
diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c
@@ -223,6 +223,50 @@ Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
 
 	PG_RETURN_POINTER(result);
 }
+
+/***********************************************************************
+ * Interpolate a point along a line 3D version
+ * --vincent.mora@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@rochgrp.com;
  ***********************************************************************/

diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
@@ -6249,6 +6249,12 @@ GRANT SELECT ON TABLE geography_columns TO public;
 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"

diff --git a/regress/core/regress_lrs.sql b/regress/core/regress_lrs.sql
@@ -72,6 +72,13 @@ select 'line_substring_12', ST_AsText(ST_LineSubstring('LINESTRING(0 0 10, 1 1 5
 select 'line_interpolate_point', ST_AsText(ST_LineInterpolatePoint('LINESTRING(0 0, 1 1)', 0));
 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
 --

diff --git a/regress/core/regress_lrs_expected b/regress/core/regress_lrs_expected
@@ -39,6 +39,8 @@ line_substring_11|POINT(0 0)
 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)