PostGIS  2.3.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 1071 of file lwlinearreferencing.c.

References compare_double(), distance3d_pt_pt(), FP_MAX, FP_MIN, lwalloc(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_as_lwline(), lwgeom_calculate_gbox(), 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  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  /* We use lwgeom_calculate_gbox() instead of lwgeom_get_gbox() */
1104  /* because we cannot afford the float rounding inaccuracy when */
1105  /* we compare the ranges for overlap below */
1106  lwgeom_calculate_gbox(g1, &gbox1);
1107  lwgeom_calculate_gbox(g2, &gbox2);
1108 
1109  /*
1110  * Find overlapping M range
1111  * WARNING: may be larger than the real one
1112  */
1113 
1114  tmin = FP_MAX(gbox1.mmin, gbox2.mmin);
1115  tmax = FP_MIN(gbox1.mmax, gbox2.mmax);
1116 
1117  if ( tmax < tmin )
1118  {
1119  LWDEBUG(1, "Inputs never exist at the same time");
1120  return -2;
1121  }
1122 
1123  // lwnotice("Min:%g, Max:%g", tmin, tmax);
1124 
1125  /*
1126  * Collect M values in common time range from inputs
1127  */
1128 
1129  mvals = lwalloc( sizeof(double) *
1130  ( l1->points->npoints + l2->points->npoints ) );
1131 
1132  /* TODO: also clip the lines ? */
1133  nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals);
1134  nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals);
1135 
1136  /* Sort values in ascending order */
1137  qsort(mvals, nmvals, sizeof(double), compare_double);
1138 
1139  /* Remove duplicated values */
1140  nmvals = uniq(mvals, nmvals);
1141 
1142  if ( nmvals < 2 )
1143  {
1144  {
1145  /* there's a single time, must be that one... */
1146  double t0 = mvals[0];
1147  POINT4D p0, p1;
1148  LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0);
1149  if ( mindist )
1150  {
1151  if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) )
1152  {
1153  lwfree(mvals);
1154  lwerror("Could not find point with M=%g on first geom", t0);
1155  return -1;
1156  }
1157  if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) )
1158  {
1159  lwfree(mvals);
1160  lwerror("Could not find point with M=%g on second geom", t0);
1161  return -1;
1162  }
1163  *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1);
1164  }
1165  lwfree(mvals);
1166  return t0;
1167  }
1168  }
1169 
1170  /*
1171  * For each consecutive pair of measures, compute time of closest point
1172  * approach and actual distance between points at that time
1173  */
1174  mintime = tmin;
1175  for (i=1; i<nmvals; ++i)
1176  {
1177  double t0 = mvals[i-1];
1178  double t1 = mvals[i];
1179  double t;
1180  POINT4D p0, p1, q0, q1;
1181  int seg;
1182  double dist2;
1183 
1184  // lwnotice("T %g-%g", t0, t1);
1185 
1186  seg = ptarray_locate_along_linear(l1->points, t0, &p0, 0);
1187  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1188  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z);
1189 
1190  seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg);
1191  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1192  // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z);
1193 
1194  seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0);
1195  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1196  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z);
1197 
1198  seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg);
1199  if ( -1 == seg ) continue; /* possible, if GBOX is approximated */
1200  // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z);
1201 
1202  t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1);
1203 
1204  /*
1205  lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g",
1206  p0.x, p0.y, p0.z,
1207  q0.x, q0.y, q0.z, t);
1208  */
1209 
1210  dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) +
1211  ( q0.y - p0.y ) * ( q0.y - p0.y ) +
1212  ( q0.z - p0.z ) * ( q0.z - p0.z );
1213  if ( dist2 < mindist2 )
1214  {
1215  mindist2 = dist2;
1216  mintime = t;
1217  // lwnotice("MINTIME: %g", mintime);
1218  }
1219  }
1220 
1221  /*
1222  * Release memory
1223  */
1224 
1225  lwfree(mvals);
1226 
1227  if ( mindist )
1228  {
1229  *mindist = sqrt(mindist2);
1230  }
1231  /*lwnotice("MINDIST: %g", sqrt(mindist2));*/
1232 
1233  return mintime;
1234 }
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)
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)
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: