Commit a714ff44 by Sandro Santilli

```Based on new lwgeom_tcpa liblwgeom function.
Includes unit and regress tests.
Includes documentation.

git-svn-id: http://svn.osgeo.org/postgis/trunk@13560 b70326c6-7e19-0410-871a-916f4a2858ee```
parent 69a0e928
 ... ... @@ -29,6 +29,7 @@ PostGIS 2.2.0 * 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) ... ...
 ... ... @@ -542,7 +542,7 @@ LINESTRING(6.1 7.1 6,7 8 9) ST_AddMeasure 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. Return a derived geometry with measure elements linearly interpolated between the start and end points. ... ... @@ -597,6 +597,84 @@ ST_GeomFromEWKT('MULTILINESTRINGM((1 0 4, 2 0 4, 4 0 4),(1 0 4, 2 0 4, 4 0 4))') ST_ClosestPointOfApproach Returns the measure at which points interpolated along two lines are closest. float8 ST_ClosestPointOfApproach geometry track1 geometry track2 Description 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. See for getting the actual points at the given measure. Availability: 2.2.0 &Z_support; Examples -- 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 See Also ,
 ... ... @@ -2,7 +2,9 @@ * * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * Copyright 2009 Paul Ramsey * * Copyright (C) 2015 Sandro Santilli * Copyright (C) 2009 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. ... ... @@ -869,6 +871,226 @@ test_lw_dist2d_ptarray_ptarrayarc(void) 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. */ ... ... @@ -887,4 +1109,5 @@ void measures_suite_setup(void) 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); }
 ... ... @@ -250,6 +250,7 @@ cu_errorreporter(const char *fmt, va_list ap) { 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 ... ...
 ... ... @@ -3,6 +3,8 @@ * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * * Copyright (C) 2009 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. * ... ... @@ -20,3 +22,11 @@ void cu_error_msg_reset(void); /* 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);
 ... ... @@ -1396,6 +1396,16 @@ extern LWCOLLECTION* lwgeom_locate_between(const LWGEOM *lwin, double from, doub */ 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. ... ...
 ... ... @@ -2,7 +2,9 @@ * * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * Copyright 2011 Paul Ramsey * * Copyright (C) 2015 Sandro Santilli * 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. ... ... @@ -11,7 +13,7 @@ #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) ... ... @@ -849,3 +851,350 @@ lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt) } 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; inpoints; ++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; inpoints; 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; ipoints, 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);*/