* New Features *
+ - #3128, ST_ClosestPointOfApproach (Sandro Santilli / Boundless)
- Canonical output for index key types
- ST_SwapOrdinates (Sandro Santilli / Boundless)
- #2918, Use GeographicLib functions for geodetics (Mike Toews)
<refnamediv>
<refname>ST_AddMeasure</refname>
- <refpurpose>Return a derived geometry with measure elements linearly interpolated between the start and end points. If the geometry has no measure dimension, one is added. If the geometry has a measure dimension, it is over-written with new values. Only LINESTRINGS and MULTILINESTRINGS are supported.</refpurpose>
+ <refpurpose>Return a derived geometry with measure elements linearly interpolated between the start and end points.</refpurpose>
</refnamediv>
<refsynopsisdiv>
</refsection>
</refentry>
-
+ <refentry id="ST_ClosestPointOfApproach">
+
+ <refnamediv>
+ <refname>ST_ClosestPointOfApproach</refname>
+ <refpurpose>
+Returns the measure at which points interpolated along two lines are closest.
+ </refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>float8 <function>ST_ClosestPointOfApproach</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>track1</parameter></paramdef>
+ <paramdef><type>geometry </type> <parameter>track2</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>
+Returns the smallest measure at which point interpolated along the given
+lines are at the smallest distance. Inputs must be LINESTRINGM geometries
+and the M value in their vertices are expected to be growing from first
+to last vertex.
+ </para>
+
+ <para>
+See <xref linkend="ST_LocateAlong" /> for getting the actual points at
+the given measure.
+ </para>
+
+ <para>Availability: 2.2.0</para>
+ <para>&Z_support;</para>
+ </refsection>
+
+
+ <refsection>
+ <title>Examples</title>
+<programlisting>
+-- Return the time in which two objects moving between 10:00 and 11:00
+-- are closest to each other and their distance at that point
+WITH inp AS ( SELECT
+ ST_AddMeasure('LINESTRING Z (0 0 0, 10 0 5)'::geometry,
+ extract(epoch from '2015-05-26 10:00'::timestamptz),
+ extract(epoch from '2015-05-26 11:00'::timestamptz)
+ ) a,
+ ST_AddMeasure('LINESTRING Z (0 2 10, 12 1 2)'::geometry,
+ extract(epoch from '2015-05-26 10:00'::timestamptz),
+ extract(epoch from '2015-05-26 11:00'::timestamptz)
+ ) b
+), cpa AS (
+ SELECT ST_ClosestPointOfApproach(a,b) m FROM inp
+), points AS (
+ SELECT ST_Force3DZ(ST_GeometryN(ST_LocateAlong(a,m),1)) pa,
+ ST_Force3DZ(ST_GeometryN(ST_LocateAlong(b,m),1)) pb
+ FROM inp, cpa
+)
+SELECT to_timestamp(m) t,
+ ST_Distance(pa,pb) distance
+FROM points, cpa;
+
+ t | distance
+-------------------------------+------------------
+ 2015-05-26 10:45:31.034483+02 | 1.96036833151395
+</programlisting>
+ </refsection>
+
+ <!-- Optionally add a "See Also" section -->
+ <refsection>
+ <title>See Also</title>
+ <para>
+<xref linkend="ST_LocateAlong" />,
+<xref linkend="ST_AddMeasure" />
+ </para>
+ </refsection>
+ </refentry>
</sect1>
*
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
- * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * Copyright (C) 2015 Sandro Santilli <strk@keybit.net>
+ * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
lwline_free(lwline1);
}
+static void
+test_lwgeom_tcpa(void)
+{
+ LWGEOM *g1, *g2;
+ double m, dist;
+
+ /* Invalid input, lack of dimensions */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING (0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
+ m = lwgeom_tcpa(g1, g2, NULL);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, -1.0);
+ CU_ASSERT_STRING_EQUAL(
+ "Both input geometries must have a measure dimension",
+ cu_error_msg
+ );
+
+ /* Invalid input, not linestrings */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("POINT M (0 0 2)", LW_PARSER_CHECK_NONE);
+ m = lwgeom_tcpa(g1, g2, NULL);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, -1.0);
+ CU_ASSERT_STRING_EQUAL(
+ "Both input geometries must be linestrings",
+ cu_error_msg
+ );
+
+ /* Invalid input, too short linestring */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(2 0 1)", LW_PARSER_CHECK_NONE);
+ dist = -77;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(dist, -77.0); /* not touched */
+ ASSERT_DOUBLE_EQUAL(m, -1.0);
+ CU_ASSERT_STRING_EQUAL(
+ "Both input lines must have at least 2 points",
+ cu_error_msg
+ );
+
+ /* Invalid input, empty linestring */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
+ m = lwgeom_tcpa(g1, g2, NULL);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, -1.0);
+ CU_ASSERT_STRING_EQUAL(
+ "Both input lines must have at least 2 points",
+ cu_error_msg
+ );
+
+ /* Invalid input, timeranges do not overlap */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
+ m = lwgeom_tcpa(g1, g2, NULL);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, -1.0);
+ CU_ASSERT_STRING_EQUAL(
+ "Inputs never exist at the same time",
+ cu_error_msg
+ );
+
+ /* One of the tracks is still, the other passes to that point */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 1.0);
+ ASSERT_DOUBLE_EQUAL(dist, 0.0);
+
+ /* One of the tracks is still, the other passes at 10 meters from point */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(-10 10 1, 10 10 5)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 3.0);
+ ASSERT_DOUBLE_EQUAL(dist, 10.0);
+
+ /* Equal tracks, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 10.0);
+ ASSERT_DOUBLE_EQUAL(dist, 0.0);
+
+ /* Reversed tracks, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(10 0 10, 0 0 20)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 15.0);
+ ASSERT_DOUBLE_EQUAL(dist, 0.0);
+
+ /* Parallel tracks, same speed, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(2 0 10, 12 0 20)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(13 0 10, 23 0 20)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 10.0);
+ ASSERT_DOUBLE_EQUAL(dist, 11.0);
+
+ /* Parallel tracks, different speed (g2 gets closer as time passes), 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 9 0 20)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 20.0);
+ ASSERT_DOUBLE_EQUAL(dist, 1.0);
+
+ /* Parallel tracks, different speed (g2 left behind as time passes), 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 6 0 20)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 10.0);
+ ASSERT_DOUBLE_EQUAL(dist, 2.0);
+
+ /* Tracks, colliding, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(5 -8 0, 5 8 10)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 5.0);
+ ASSERT_DOUBLE_EQUAL(dist, 0.0);
+
+ /* Tracks crossing, NOT colliding, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(8 -5 0, 8 5 10)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 6.5);
+ ASSERT_DOUBLE_EQUAL(rint(dist*100), 212.0);
+
+ /* Same origin, different direction, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, -100 0 10)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 1.0);
+ ASSERT_DOUBLE_EQUAL(dist, 0.0);
+
+ /* Same ending, different direction, 2d */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(10 0 1, 0 0 10)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(0 -100 1, 0 0 10)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 10.0);
+ ASSERT_DOUBLE_EQUAL(dist, 0.0);
+
+ /* Converging tracks, 3d */
+
+ g1 = lwgeom_from_wkt("LINESTRING ZM(0 0 0 10, 10 0 0 20)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING ZM(0 0 8 10, 10 0 5 20)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 20.0);
+ ASSERT_DOUBLE_EQUAL(dist, 5.0);
+
+
+ /* G1 stops at t=1 until t=4 to let G2 pass by, then continues */
+ /* G2 passes at 1 meter from G1 t=3 */
+
+ g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 0 1 1, 0 1 4, 0 10 13)", LW_PARSER_CHECK_NONE);
+ g2 = lwgeom_from_wkt("LINESTRING M(-10 2 0, 0 2 3, 12 2 13)", LW_PARSER_CHECK_NONE);
+ dist = -1;
+ m = lwgeom_tcpa(g1, g2, &dist);
+ lwgeom_free(g1);
+ lwgeom_free(g2);
+ ASSERT_DOUBLE_EQUAL(m, 3.0);
+ ASSERT_DOUBLE_EQUAL(dist, 1.0);
+
+}
+
/*
** Used by test harness to register the tests in this file.
*/
PG_ADD_TEST(suite, test_lw_arc_length);
PG_ADD_TEST(suite, test_lw_dist2d_pt_ptarrayarc);
PG_ADD_TEST(suite, test_lw_dist2d_ptarray_ptarrayarc);
+ PG_ADD_TEST(suite, test_lwgeom_tcpa);
}
{
vsnprintf (cu_error_msg, MAX_CUNIT_MSG_LENGTH, fmt, ap);
cu_error_msg[MAX_CUNIT_MSG_LENGTH]='\0';
+ /* fprintf(stderr, "ERROR: %s\n", cu_error_msg); */
}
void
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
*
+ * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
*
/* Our internal callback to register Suites with the main tester */
typedef void (*PG_SuiteSetup)(void);
+
+#define ASSERT_DOUBLE_EQUAL(o,e) do { \
+ if ( o != e ) \
+ fprintf(stderr, "[%s:%d]\n Expected: %g\n Obtained: %g\n", __FILE__, __LINE__, (e), (o)); \
+ CU_ASSERT_EQUAL(o,e); \
+} while (0);
+
+
*/
extern double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt);
+/**
+* Find the time of closest point of approach
+*
+* @param mindist if not null will be set to the minimum distance between
+* the trajectories at the closest point of approach.
+*
+* @return the time value in which the minimum distance was reached
+*/
+extern double lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist);
+
/*
* Ensure every segment is at most 'dist' long.
* Returned LWGEOM might is unchanged if a POINT.
*
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
- * Copyright 2011 Paul Ramsey
+ *
+ * Copyright (C) 2015 Sandro Santilli <strk@keybit.net>
+ * Copyright (C) 2011 Paul Ramsey
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
-
+#include "measures3d.h"
static int
segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn)
}
return ret;
}
+
+/*
+ * Time of closest point of approach
+ *
+ * Given two vectors (p1-p2 and q1-q2) and
+ * a time range (t1-t2) return the time in which
+ * a point p is closest to a point q on their
+ * respective vectors, and the actual points
+ *
+ * Here we use algorithm from softsurfer.com
+ * that can be found here
+ * http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
+ *
+ * @param p0 start of first segment, will be set to actual
+ * closest point of approach on segment.
+ * @param p1 end of first segment
+ * @param q0 start of second segment, will be set to actual
+ * closest point of approach on segment.
+ * @param q1 end of second segment
+ * @param t0 start of travel time
+ * @param t1 end of travel time
+ *
+ * @return time of closest point of approach
+ *
+ */
+static double
+segments_tcpa(POINT4D* p0, const POINT4D* p1,
+ POINT4D* q0, const POINT4D* q1,
+ double t0, double t1)
+{
+ POINT3DZ pv; /* velocity of p, aka u */
+ POINT3DZ qv; /* velocity of q, aka v */
+ POINT3DZ dv; /* velocity difference */
+ POINT3DZ w0; /* vector between first points */
+
+/*
+ lwnotice("FROM %g,%g,%g,%g -- %g,%g,%g,%g",
+ p0->x, p0->y, p0->z, p0->m,
+ p1->x, p1->y, p1->z, p1->m);
+ lwnotice(" TO %g,%g,%g,%g -- %g,%g,%g,%g",
+ q0->x, q0->y, q0->z, q0->m,
+ q1->x, q1->y, q1->z, q1->m);
+*/
+
+ /* PV aka U */
+ pv.x = ( p1->x - p0->x );
+ pv.y = ( p1->y - p0->y );
+ pv.z = ( p1->z - p0->z );
+ /*lwnotice("PV: %g, %g, %g", pv.x, pv.y, pv.z);*/
+
+ /* QV aka V */
+ qv.x = ( q1->x - q0->x );
+ qv.y = ( q1->y - q0->y );
+ qv.z = ( q1->z - q0->z );
+ /*lwnotice("QV: %g, %g, %g", qv.x, qv.y, qv.z);*/
+
+ dv.x = pv.x - qv.x;
+ dv.y = pv.y - qv.y;
+ dv.z = pv.z - qv.z;
+ /*lwnotice("DV: %g, %g, %g", dv.x, dv.y, dv.z);*/
+
+ double dv2 = DOT(dv,dv);
+ /*lwnotice("DOT: %g", dv2);*/
+
+ if ( dv2 == 0.0 )
+ {
+ /* Distance is the same at any time, we pick the earliest */
+ return t0;
+ }
+
+ /* Distance at any given time, with t0 */
+ w0.x = ( p0->x - q0->x );
+ w0.y = ( p0->y - q0->y );
+ w0.z = ( p0->z - q0->z );
+
+ /*lwnotice("W0: %g, %g, %g", w0.x, w0.y, w0.z);*/
+
+ /* Check that at distance dt w0 is distance */
+
+ /* This is the fraction of measure difference */
+ double t = -DOT(w0,dv) / dv2;
+ /*lwnotice("CLOSEST TIME (fraction): %g", t);*/
+
+ if ( t > 1.0 ) {
+ /* Getting closer as we move to the end */
+ /*lwnotice("Converging");*/
+ t = 1;
+ } else if ( t < 0.0 ) {
+ /*lwnotice("Diverging");*/
+ t = 0;
+ }
+
+ /* Interpolate the actual points now */
+
+ p0->x += pv.x * t;
+ p0->y += pv.y * t;
+ p0->z += pv.z * t;
+
+ q0->x += qv.x * t;
+ q0->y += qv.y * t;
+ q0->z += qv.z * t;
+
+ t = t0 + (t1 - t0) * t;
+ /*lwnotice("CLOSEST TIME (real): %g", t);*/
+
+ return t;
+}
+
+static int
+ptarray_collect_mvals(const POINTARRAY *pa, double tmin, double tmax, double *mvals)
+{
+ POINT4D pbuf;
+ int i, n=0;
+ for (i=0; i<pa->npoints; ++i)
+ {
+ getPoint4d_p(pa, i, &pbuf); /* could be optimized */
+ if ( pbuf.m >= tmin && pbuf.m <= tmax )
+ mvals[n++] = pbuf.m;
+ }
+ return n;
+}
+
+static int
+compare_double(const void *a, const void *b)
+{
+ double *dpa = (double *)a;
+ double *dpb = (double *)b;
+ return *dpb < *dpa;
+}
+
+/* Return number of elements in unique array */
+static int
+uniq(double *vals, int nvals)
+{
+ int i, last=0;
+ for (i=1; i<nvals; ++i)
+ {
+ /*lwnotice("(I%d):%g", i, vals[i]);*/
+ if ( vals[i] != vals[last] )
+ {
+ vals[++last] = vals[i];
+ /*lwnotice("(O%d):%g", last, vals[last]);*/
+ }
+ }
+ return last+1;
+}
+
+/*
+ * Find point at a given measure
+ *
+ * The function assumes measures are linear so that always a single point
+ * is returned for a single measure.
+ *
+ * @param pa the point array to perform search on
+ * @param m the measure to search for
+ * @param p the point to write result into
+ * @param from the segment number to start from
+ *
+ * @return the segment number the point was found into
+ * or -1 if given measure was out of the known range.
+ */
+static int
+ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from)
+{
+ int i = from;
+ POINT4D p1, p2;
+
+ /* Walk through each segment in the point array */
+ getPoint4d_p(pa, i, &p1);
+ for ( i = 1; i < pa->npoints; i++ )
+ {
+ getPoint4d_p(pa, i, &p2);
+
+ if ( segment_locate_along(&p1, &p2, m, 0, p) == LW_TRUE )
+ return i-1; /* found */
+
+ p1 = p2;
+ }
+
+ return -1; /* not found */
+}
+
+double
+lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist)
+{
+ LWLINE *l1, *l2;
+ int i;
+ const GBOX *gbox1, *gbox2;
+ double tmin, tmax;
+ double *mvals;
+ int nmvals = 0;
+ double mintime;
+ double mindist2; /* minimum distance, squared */
+
+ if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) ) {
+ lwerror("Both input geometries must have a measure dimension");
+ return -1;
+ }
+
+ l1 = lwgeom_as_lwline(g1);
+ l2 = lwgeom_as_lwline(g2);
+
+ if ( ! l1 || ! l2 ) {
+ lwerror("Both input geometries must be linestrings");
+ return -1;
+ }
+
+ if ( l1->points->npoints < 2 || l2->points->npoints < 2 ) {
+ lwerror("Both input lines must have at least 2 points");
+ return -1;
+ }
+
+ gbox1 = lwgeom_get_bbox(g1);
+ gbox2 = lwgeom_get_bbox(g2);
+
+ assert(gbox1); /* or the npoints check above would have failed */
+ assert(gbox2); /* or the npoints check above would have failed */
+
+ /*
+ * Find overlapping M range
+ */
+
+ tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
+ tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
+
+ if ( tmax == tmin ) /* both exists only at a given time */
+ {
+ /*lwnotice("Inputs only exist both at a single time");*/
+ if ( mindist )
+ {
+ *mindist = lwgeom_mindistance3d(g1, g2);
+ }
+ return tmax;
+ }
+
+ if ( tmax < tmin ) {
+ lwerror("Inputs never exist at the same time");
+ return -1;
+ }
+
+ /* lwnotice("Min:%g, Max:%g", tmin, tmax); */
+
+ /*
+ * Collect M values in common time range from inputs
+ */
+
+ mvals = lwalloc( sizeof(double*) *
+ ( l1->points->npoints + l2->points->npoints ) );
+
+ /* TODO: also clip the lines ? */
+ nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
+ nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, &(mvals[nmvals]));
+
+ /* Sort values in ascending order */
+ qsort(mvals, nmvals, sizeof(double), compare_double);
+
+ /* Remove duplicated values */
+ nmvals = uniq(mvals, nmvals);
+
+ if ( nmvals < 2 )
+ return mvals[0]; /* there's a single time, must be that one... */
+
+ /*
+ * For each consecutive pair of measures, compute time of closest point
+ * approach and actual distance between points at that time
+ */
+ mintime = tmin;
+ for (i=1; i<nmvals; ++i)
+ {
+ double t0 = mvals[i-1];
+ double t1 = mvals[i];
+ double t;
+ POINT4D p0, p1, q0, q1;
+ int seg;
+ double dist2;
+
+ /*lwnotice("T %g-%g", t0, t1);*/
+
+ seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
+ if ( -1 == seg )
+ {
+ lwfree(mvals);
+ lwerror("Non-linear measures in first geometry");
+ return -1;
+ }
+ /*lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);*/
+
+ seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
+ if ( -1 == seg )
+ {
+ lwfree(mvals);
+ lwerror("Non-linear measures in first geometry");
+ return -1;
+ }
+ /*lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);*/
+
+ seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
+ if ( -1 == seg )
+ {
+ lwfree(mvals);
+ lwerror("Non-linear measures in second geometry");
+ return -1;
+ }
+ /*lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);*/
+
+ seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
+ if ( -1 == seg )
+ {
+ lwfree(mvals);
+ lwerror("Non-linear measures in second geometry");
+ return -1;
+ }
+ /*lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);*/
+
+ t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
+
+/*
+ lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
+ p0.x, p0.y, p0.z,
+ q0.x, q0.y, q0.z, t);
+*/
+
+ dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
+ ( q0.y - p0.y ) * ( q0.y - p0.y ) +
+ ( q0.z - p0.z ) * ( q0.z - p0.z );
+ if ( i == 1 || dist2 < mindist2 )
+ {
+ mindist2 = dist2;
+ mintime = t;
+ /* lwnotice("MINTIME: %g", mintime); */
+ }
+
+ }
+
+ /*
+ * Release memory
+ */
+
+ lwfree(mvals);
+
+ if ( mindist ) {
+ *mindist = sqrt(mindist2);
+ }
+ /*lwnotice("MINDIST: %g", sqrt(mindist2));*/
+
+ return mintime;
+}
PG_RETURN_POINTER(gout);
}
+/***********************************************************************
+ *
+ * Spatio-Temporal functions
+ *
+ ************************************************************************/
+
+/*
+ * Return the measure at which interpolated points on the two
+ * input lines are at the smallest distance.
+ */
+Datum ST_ClosestPointOfApproach(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_ClosestPointOfApproach);
+Datum ST_ClosestPointOfApproach(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gs0 = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *gs1 = PG_GETARG_GSERIALIZED_P(1);
+ /* All checks already performed by liblwgeom, not worth checking again */
+ LWGEOM *g0 = lwgeom_from_gserialized(gs0);
+ LWGEOM *g1 = lwgeom_from_gserialized(gs1);
+ double m = lwgeom_tcpa(g0, g1, NULL);
+ lwgeom_free(g0);
+ lwgeom_free(g1);
+ PG_FREE_IF_COPY(gs0, 0);
+ PG_FREE_IF_COPY(gs1, 1);
+ PG_RETURN_FLOAT8(m);
+}
RETURNS geometry
AS 'MODULE_PATHNAME', 'ST_AddMeasure'
LANGUAGE 'c' IMMUTABLE STRICT;
+
+---------------------------------------------------------------
+-- TEMPORAL
+---------------------------------------------------------------
+
+-- Availability: 2.2.0
+CREATE OR REPLACE FUNCTION ST_ClosestPointOfApproach(geometry, geometry)
+ RETURNS float8
+ AS 'MODULE_PATHNAME', 'ST_ClosestPointOfApproach'
+ LANGUAGE 'c' IMMUTABLE STRICT;
---------------------------------------------------------------
-- GEOS
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_PointOfClosestApproach
+
+-- Converging
+select 'pca1', ST_ClosestPointOfApproach(
+ 'LINESTRINGZM(0 0 0 0, 10 10 10 10)',
+ 'LINESTRINGZM(0 0 0 1, 10 10 10 10)');
+-- Following
+select 'pca2', ST_ClosestPointOfApproach(
+ 'LINESTRINGZM(0 0 0 0, 10 10 10 10)',
+ 'LINESTRINGZM(0 0 0 5, 10 10 10 15)');
+-- Crossing
+select 'pca3', ST_ClosestPointOfApproach(
+ 'LINESTRINGZM(0 0 0 0, 0 0 0 10)',
+ 'LINESTRINGZM(-30 0 5 4, 10 0 5 6)');
+
+-- TODO: test ST_AddMeasures
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)
+pca1|10
+pca2|5
+pca3|5.5
projects / postgis / commitdiff