PostGIS  2.2.8dev-r@@SVN_REVISION@@

◆ lwgeom_tcpa()

double lwgeom_tcpa ( const LWGEOM g1,
const LWGEOM g2,
double *  mindist 
)

Find the time of closest point of approach.

Parameters
mindistif not null will be set to the minimum distance between the trajectories at the closest point of approach.
Returns
the time value in which the minimum distance was reached, -1 if inputs are invalid (lwerror is called in that case), -2 if the trajectories do not share any point in time.

Definition at line 1058 of file lwlinearreferencing.c.

References compare_double(), distance3d_pt_pt(), FP_MAX, FP_MIN, lwalloc(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_as_lwline(), lwgeom_get_bbox(), lwgeom_has_m(), 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_ClosestPointOfApproach(), ST_DistanceCPA(), and test_lwgeom_tcpa().

1059 {
1060  LWLINE *l1, *l2;
1061  int i;
1062  const GBOX *gbox1, *gbox2;
1063  double tmin, tmax;
1064  double *mvals;
1065  int nmvals = 0;
1066  double mintime;
1067  double mindist2 = FLT_MAX; /* minimum distance, squared */
1068 
1069  if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) )
1070  {
1071  lwerror("Both input geometries must have a measure dimension");
1072  return -1;
1073  }
1074 
1075  l1 = lwgeom_as_lwline(g1);
1076  l2 = lwgeom_as_lwline(g2);
1077 
1078  if ( ! l1 || ! l2 )
1079  {
1080  lwerror("Both input geometries must be linestrings");
1081  return -1;
1082  }
1083 
1084  if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
1085  {
1086  lwerror("Both input lines must have at least 2 points");
1087  return -1;
1088  }
1089 
1090  /* WARNING: these ranges may be wider than real ones */
1091  gbox1 = lwgeom_get_bbox(g1);
1092  gbox2 = lwgeom_get_bbox(g2);
1093 
1094  assert(gbox1); /* or the npoints check above would have failed */
1095  assert(gbox2); /* or the npoints check above would have failed */
1096 
1097  /*
1098  * Find overlapping M range
1099  * WARNING: may be larger than the real one
1100  */
1101 
1102  tmin = FP_MAX(gbox1->mmin, gbox2->mmin);
1103  tmax = FP_MIN(gbox1->mmax, gbox2->mmax);
1104 
1105  if ( tmax < tmin )
1106  {
1107  LWDEBUG(1, "Inputs never exist at the same time");
1108  return -2;
1109  }
1110 
1111  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1112 
1113  /*
1114  * Collect M values in common time range from inputs
1115  */
1116 
1117  mvals = lwalloc( sizeof(double) *
1118  ( l1->points->npoints + l2->points->npoints ) );
1119 
1120  /* TODO: also clip the lines ? */
1121  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1122  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1123 
1124  /* Sort values in ascending order */
1125  qsort(mvals, nmvals, sizeof(double), compare_double);
1126 
1127  /* Remove duplicated values */
1128  nmvals = uniq(mvals, nmvals);
1129 
1130  if ( nmvals < 2 )
1131  {
1132  {
1133  /* there's a single time, must be that one... */
1134  double t0 = mvals[0];
1135  POINT4D p0, p1;
1136  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1137  if ( mindist )
1138  {
1139  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1140  {
1141  lwfree(mvals);
1142  lwerror("Could not find point with M=%g on first geom", t0);
1143  return -1;
1144  }
1145  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1146  {
1147  lwfree(mvals);
1148  lwerror("Could not find point with M=%g on second geom", t0);
1149  return -1;
1150  }
1151  *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1);
1152  }
1153  lwfree(mvals);
1154  return t0;
1155  }
1156  }
1157 
1158  /*
1159  * For each consecutive pair of measures, compute time of closest point
1160  * approach and actual distance between points at that time
1161  */
1162  mintime = tmin;
1163  for (i=1; i<nmvals; ++i)
1164  {
1165  double t0 = mvals[i-1];
1166  double t1 = mvals[i];
1167  double t;
1168  POINT4D p0, p1, q0, q1;
1169  int seg;
1170  double dist2;
1171 
1172  // lwnotice("T %g-%g", t0, t1);
1173 
1174  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1175  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1176  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1177 
1178  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1179  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1180  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1181 
1182  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1183  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1184  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1185 
1186  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1187  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1188  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1189 
1190  t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1191 
1192  /*
1193  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1194  p0.x, p0.y, p0.z,
1195  q0.x, q0.y, q0.z, t);
1196  */
1197 
1198  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1199  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1200  ( q0.z - p0.z ) * ( q0.z - p0.z );
1201  if ( dist2 < mindist2 )
1202  {
1203  mindist2 = dist2;
1204  mintime = t;
1205  // lwnotice("MINTIME: %g", mintime);
1206  }
1207  }
1208 
1209  /*
1210  * Release memory
1211  */
1212 
1213  lwfree(mvals);
1214 
1215  if ( mindist )
1216  {
1217  *mindist = sqrt(mindist2);
1218  }
1219  /*lwnotice("MINDIST: %g", sqrt(mindist2));*/
1220 
1221  return mintime;
1222 }
double x
Definition: liblwgeom.h:336
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)
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: