PostGIS  2.3.7dev-r@@SVN_REVISION@@
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 1072 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().

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