PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ lwgeom_cpa_within()

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 1237 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_calculate_gbox(), 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().

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