PostGIS  2.5.0beta1dev-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 1243 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().

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