PostGIS  2.5.0dev-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 1238 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: