[postgis-devel] Interpolating a point along a line

[Joe Sunday]

On Wed, Jan 21, 2004 at 09:00:39AM -0800, David Blasby wrote:

> This is a good function - but I think we should make 2 minor changes to 
> strk to apply the patch and make sure its working.

Corrected patch applied. 

--Joe

-- 
Joe Sunday <sunday at csh.rit.edu>  http://www.csh.rit.edu/~sunday/
Computer Science House, Rochester Inst. Of Technology
-------------- next part --------------
Index: postgis.h
===================================================================
RCS file: /home/cvs/postgis/postgis/postgis.h,v
retrieving revision 1.39
diff -u -5 -p -r1.39 postgis.h
--- postgis.h	28 Nov 2003 11:06:49 -0000	1.39
+++ postgis.h	21 Jan 2004 18:04:21 -0000
@@ -631,10 +631,11 @@ Datum geometry_from_text_mpoint(PG_FUNCT
 Datum geometry_from_text_line(PG_FUNCTION_ARGS);
 Datum geometry_from_text_mline(PG_FUNCTION_ARGS);
 Datum geometry_from_text_gc(PG_FUNCTION_ARGS);
 Datum isempty(PG_FUNCTION_ARGS);
 Datum simplify(PG_FUNCTION_ARGS);
+Datum line_interpolate_point(PG_FUNCTION_ARGS);
 
 /*--------------------------------------------------------------------
  * Useful floating point utilities and constants.
  * from postgres geo_decls.c
  * EPSILON modified to be more "double" friendly
Index: postgis_algo.c
===================================================================
RCS file: /home/cvs/postgis/postgis/postgis_algo.c,v
retrieving revision 1.1
diff -u -5 -p -r1.1 postgis_algo.c
--- postgis_algo.c	27 Oct 2003 10:21:15 -0000	1.1
+++ postgis_algo.c	21 Jan 2004 18:04:21 -0000
@@ -405,5 +405,101 @@ Datum simplify(PG_FUNCTION_ARGS)
 }
 
 /***********************************************************************
  * --strk at keybit.net;
  ***********************************************************************/
+
+/***********************************************************************
+ * Interpolate a point along a line, useful for Geocoding applications
+ * SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
+ * returns POINT( 1 1 )
+ *
+ * -- jsunday at rochgrp.com;
+ ***********************************************************************/
+PG_FUNCTION_INFO_V1(line_interpolate_point);
+Datum line_interpolate_point(PG_FUNCTION_ARGS)
+{
+	GEOMETRY *geom = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	double distance = PG_GETARG_FLOAT8(1);
+
+	int32 *offsets1;
+	LINE3D *line;
+	POINT3D point;
+	int nsegs, i;
+	double length, slength, tlength;
+
+	if( distance < 0 || distance > 1 ) {
+		elog(ERROR,"line_interpolate_point: 2nd arg isnt within [0,1]");
+		PG_RETURN_NULL();
+	}
+
+	if( geom->type != LINETYPE ) {
+		elog(ERROR,"line_interpolate_point: 1st arg isnt a line");
+		PG_RETURN_NULL();
+	}
+
+	offsets1 = (int32 *) ( ((char *) &(geom->objType[0] )) +
+			sizeof(int32) * geom->nobjs );
+	line = (LINE3D*) ( (char *)geom + offsets1[0] );
+
+	/* If distance is one of the two extremes, return the point on that
+	 * end rather than doing any expensive computations
+	 */
+	if( distance == 0.0 ) {
+		PG_RETURN_POINTER(
+				make_oneobj_geometry(sizeof(POINT3D),
+					/*(char*)&(line->points[0]), POINTTYPE, */
+					(char*)&(line->points[0]), POINTTYPE,
+					geom->is3d, geom->SRID, geom->scale, geom->offsetX, geom->offsetY
+					)
+				);
+	}
+
+	if( distance == 1.0 ) {
+		PG_RETURN_POINTER(
+				make_oneobj_geometry(sizeof(POINT3D),
+					(char*)&(line->points[line->npoints-1]), POINTTYPE, 
+					geom->is3d, geom->SRID, geom->scale, geom->offsetX, geom->offsetY
+					)
+				);
+	}
+
+	/* Interpolate a point on the line */
+	nsegs = line->npoints - 1;
+	length = line_length2d( line );
+	tlength = 0;
+	for( i = 0; i < nsegs; i++ ) {
+		POINT3D *p1, *p2;
+		p1 = &(line->points[i]);
+		p2 = &(line->points[i+1]);
+		/* Find the relative length of this segment */
+		slength = distance_pt_pt( p1, p2 )/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;
+			point.x = (p1->x) + ((p2->x - p1->x) * dseg);
+			point.y = (p1->y) + ((p2->y - p1->y) * dseg);
+			point.z = 0;
+			PG_RETURN_POINTER(
+					make_oneobj_geometry(sizeof(POINT3D),
+						(char*)&point, POINTTYPE, 
+						geom->is3d, geom->SRID, geom->scale, geom->offsetX, geom->offsetY
+						)
+					);
+		}
+		tlength += slength;
+	}
+	/* Return the last point on the line. This shouldn't happen, but
+	 * could if there's some floating point rounding errors. */
+	PG_RETURN_POINTER(
+			make_oneobj_geometry(sizeof(POINT3D),
+				(char*)&(line->points[line->npoints-1]), POINTTYPE, 
+				geom->is3d, geom->SRID, geom->scale, geom->offsetX, geom->offsetY
+				)
+			);
+}
+/***********************************************************************
+ * --jsunday at rochgrp.com;
+ ***********************************************************************/
Index: postgis_sql_common.sql.in
===================================================================
RCS file: /home/cvs/postgis/postgis/postgis_sql_common.sql.in,v
retrieving revision 1.26
diff -u -5 -p -r1.26 postgis_sql_common.sql.in
--- postgis_sql_common.sql.in	23 Dec 2003 09:00:12 -0000	1.26
+++ postgis_sql_common.sql.in	21 Jan 2004 18:04:24 -0000
@@ -1243,5 +1243,11 @@ CREATE FUNCTION geosnoop(geometry)
 
 CREATE FUNCTION simplify(geometry, float8)
    RETURNS geometry
    AS '@MODULE_FILENAME@'
    LANGUAGE 'C' WITH (isstrict);
+
+CREATE FUNCTION line_interpolate_point(geometry, float8)
+   RETURNS geometry
+	 AS '@MODULE_FILENAME@'
+	 LANGUAGE 'C' WITH (isstrict);
+