diff --git a/lwgeom/liblwgeom.h b/lwgeom/liblwgeom.h
index 3ffc0ad80..c5ced5b4b 100644
--- a/lwgeom/liblwgeom.h
+++ b/lwgeom/liblwgeom.h
@@ -934,6 +934,8 @@ extern int lwpoly_compute_box2d_p(LWPOLY *poly, BOX2DFLOAT4 *box);
 extern int lwcollection_compute_box2d_p(LWCOLLECTION *col, BOX2DFLOAT4 *box);
 extern BOX2DFLOAT4 *lwgeom_compute_box2d(LWGEOM *lwgeom);
 
+extern void interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F);
+
 /* return alloced memory */
 extern BOX2DFLOAT4 *box2d_union(BOX2DFLOAT4 *b1, BOX2DFLOAT4 *b2);
 
diff --git a/lwgeom/lwgeom_api.c b/lwgeom/lwgeom_api.c
index 2d638fd58..f11d36ba4 100644
--- a/lwgeom/lwgeom_api.c
+++ b/lwgeom/lwgeom_api.c
@@ -475,12 +475,11 @@ getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *op)
 	int zmflag;
 
 #if PARANOIA_LEVEL > 0
-	if ( ! pa ) return 0;
+	if ( ! pa ) lwerror("getPoint4d_p: NULL pointarray");
 
 	if ( (n<0) || (n>=pa->npoints))
 	{
 		lwerror("getPoint4d_p: point offset out of range");
-		return 0; /*error */
 	}
 #endif
 
@@ -2135,3 +2134,32 @@ parse_lwgeom_wkt(char *wkt_input)
 	return serialized_form;
 
 }
+
+/*
+ * Find interpolation point I
+ * between point A and point B 
+ * so that the len(AI) == len(AB)*F
+ * and I falls on AB segment.
+ *
+ * Example:
+ *
+ *   F=0.5  :    A----I----B 
+ *   F=1    :    A---------B==I 
+ *   F=0    : A==I---------B
+ *   F=.2   :    A-I-------B
+ */
+void
+interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F)
+{
+#if PARANOIA_LEVEL > 0
+	double absF=fabs(F);
+	if ( absF < 0 || absF > 1 )
+	{
+		lwerror("interpolate_point4d: invalid F (%g)", F);
+	}
+#endif
+	I->x=A->x+((B->x-A->x)*F);
+	I->y=A->y+((B->y-A->y)*F);
+	I->z=A->z+((B->z-A->z)*F);
+	I->m=A->m+((B->m-A->m)*F);
+}