PostGIS  2.2.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 1225 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: