PostGIS  2.5.0dev-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 1071 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: