PostGIS  2.5.0beta2dev-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 1077 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().

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