PostGIS  2.3.7dev-r@@SVN_REVISION@@
int lwgeom_cpa_within ( const LWGEOM g1,
const LWGEOM g2,
double  maxdist 
)

Is the closest point of approach within a distance ?

Returns
LW_TRUE or LW_FALSE

Definition at line 1239 of file lwlinearreferencing.c.

References compare_double(), distance3d_pt_pt(), FP_MAX, FP_MIN, LW_FALSE, LW_TRUE, lwalloc(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_as_lwline(), lwgeom_get_bbox(), lwgeom_has_m(), lwnotice(), GBOX::mmax, GBOX::mmin, POINTARRAY::npoints, LWLINE::points, ptarray_collect_mvals(), ptarray_locate_along_linear(), segments_tcpa(), uniq(), POINT4D::x, POINT4D::y, and POINT4D::z.

Referenced by ST_CPAWithin(), and test_lwgeom_tcpa().

1240 {
1241  LWLINE *l1, *l2;
1242  int i;
1243  const GBOX *gbox1, *gbox2;
1244  double tmin, tmax;
1245  double *mvals;
1246  int nmvals = 0;
1247  double maxdist2 = maxdist * maxdist;
1248  int within = LW_FALSE;
1249 
1250  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1251  {
1252  lwerror("Both input geometries must have a measure dimension");
1253  return LW_FALSE;
1254  }
1255 
1256  l1 = lwgeom_as_lwline(g1);
1257  l2 = lwgeom_as_lwline(g2);
1258 
1259  if ( ! l1 || ! l2 )
1260  {
1261  lwerror("Both input geometries must be linestrings");
1262  return LW_FALSE;
1263  }
1264 
1265  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1266  {
1267  /* TODO: return distance between these two points */
1268  lwerror("Both input lines must have at least 2 points");
1269  return LW_FALSE;
1270  }
1271 
1272  /* WARNING: these ranges may be wider than real ones */
1273  gbox1 = lwgeom_get_bbox(g1);
1274  gbox2 = lwgeom_get_bbox(g2);
1275 
1276  assert(gbox1); /* or the npoints check above would have failed */
1277  assert(gbox2); /* or the npoints check above would have failed */
1278 
1279  /*
1280  * Find overlapping M range
1281  * WARNING: may be larger than the real one
1282  */
1283 
1284  tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
1285  tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
1286 
1287  if ( tmax < tmin )
1288  {
1289  LWDEBUG(1, "Inputs never exist at the same time");
1290  return LW_FALSE;
1291  }
1292 
1293  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1294 
1295  /*
1296  * Collect M values in common time range from inputs
1297  */
1298 
1299  mvals = lwalloc( sizeof(double) *
1300  ( l1->points->npoints + l2->points->npoints ) );
1301 
1302  /* TODO: also clip the lines ? */
1303  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1304  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1305 
1306  /* Sort values in ascending order */
1307  qsort(mvals, nmvals, sizeof(double), compare_double);
1308 
1309  /* Remove duplicated values */
1310  nmvals = uniq(mvals, nmvals);
1311 
1312  if ( nmvals < 2 )
1313  {
1314  /* there's a single time, must be that one... */
1315  double t0 = mvals[0];
1316  POINT4D p0, p1;
1317  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1318  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1319  {
1320  lwnotice("Could not find point with M=%g on first geom", t0);
1321  return LW_FALSE;
1322  }
1323  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1324  {
1325  lwnotice("Could not find point with M=%g on second geom", t0);
1326  return LW_FALSE;
1327  }
1328  if ( distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1) <= maxdist )
1329  within = LW_TRUE;
1330  lwfree(mvals);
1331  return within;
1332  }
1333 
1334  /*
1335  * For each consecutive pair of measures, compute time of closest point
1336  * approach and actual distance between points at that time
1337  */
1338  for (i=1; i<nmvals; ++i)
1339  {
1340  double t0 = mvals[i-1];
1341  double t1 = mvals[i];
1342 #if POSTGIS_DEBUG_LEVEL >= 1
1343  double t;
1344 #endif
1345  POINT4D p0, p1, q0, q1;
1346  int seg;
1347  double dist2;
1348 
1349  // lwnotice("T %g-%g", t0, t1);
1350 
1351  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1352  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1353  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1354 
1355  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1356  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1357  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1358 
1359  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1360  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1361  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1362 
1363  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1364  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1365  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1366 
1367 #if POSTGIS_DEBUG_LEVEL >= 1
1368  t =
1369 #endif
1370  segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1371 
1372  /*
1373  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1374  p0.x, p0.y, p0.z,
1375  q0.x, q0.y, q0.z, t);
1376  */
1377 
1378  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1379  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1380  ( q0.z - p0.z ) * ( q0.z - p0.z );
1381  if ( dist2 <= maxdist2 )
1382  {
1383  LWDEBUGF(1, "Within distance %g at time %g, breaking", sqrt(dist2), t);
1384  within = LW_TRUE;
1385  break;
1386  }
1387  }
1388 
1389  /*
1390  * Release memory
1391  */
1392 
1393  lwfree(mvals);
1394 
1395  return within;
1396 }
double x
Definition: liblwgeom.h:351
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:89
void lwfree(void *mem)
Definition: lwutil.c:242
int npoints
Definition: liblwgeom.h:370
static int uniq(double *vals, int nvals)
static int compare_double(const void *pa, const void *pb)
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define FP_MIN(A, B)
static int ptarray_collect_mvals(const POINTARRAY *pa, double tmin, double tmax, double *mvals)
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
const GBOX * lwgeom_get_bbox(const LWGEOM *lwgeom)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:653
static int ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from)
double z
Definition: liblwgeom.h:351
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:102
double mmin
Definition: liblwgeom.h:297
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:795
double mmax
Definition: liblwgeom.h:298
void * lwalloc(size_t size)
Definition: lwutil.c:227
double y
Definition: liblwgeom.h:351
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:856
static double segments_tcpa(POINT4D *p0, const POINT4D *p1, POINT4D *q0, const POINT4D *q1, double t0, double t1)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
#define FP_MAX(A, B)
POINTARRAY * points
Definition: liblwgeom.h:421

Here is the call graph for this function:

Here is the caller graph for this function: