PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ edge_intersection()

int edge_intersection ( const GEOGRAPHIC_EDGE e1,
const GEOGRAPHIC_EDGE e2,
GEOGRAPHIC_POINT g 
)

Returns true if an intersection can be calculated, and places it in *g.

Returns false otherwise.

Definition at line 1127 of file lwgeodetic.c.

1128 {
1129  POINT3D ea, eb, v;
1130  LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", e1->start.lat, e1->start.lon, e1->end.lat, e1->end.lon);
1131  LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", e2->start.lat, e2->start.lon, e2->end.lat, e2->end.lon);
1132 
1133  LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e1->start.lon), rad2deg(e1->start.lat), rad2deg(e1->end.lon), rad2deg(e1->end.lat));
1134  LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e2->start.lon), rad2deg(e2->start.lat), rad2deg(e2->end.lon), rad2deg(e2->end.lat));
1135 
1136  if ( geographic_point_equals(&(e1->start), &(e2->start)) )
1137  {
1138  *g = e1->start;
1139  return LW_TRUE;
1140  }
1141  if ( geographic_point_equals(&(e1->end), &(e2->end)) )
1142  {
1143  *g = e1->end;
1144  return LW_TRUE;
1145  }
1146  if ( geographic_point_equals(&(e1->end), &(e2->start)) )
1147  {
1148  *g = e1->end;
1149  return LW_TRUE;
1150  }
1151  if ( geographic_point_equals(&(e1->start), &(e2->end)) )
1152  {
1153  *g = e1->start;
1154  return LW_TRUE;
1155  }
1156 
1157  robust_cross_product(&(e1->start), &(e1->end), &ea);
1158  normalize(&ea);
1159  robust_cross_product(&(e2->start), &(e2->end), &eb);
1160  normalize(&eb);
1161  LWDEBUGF(4, "e1 cross product == POINT(%.12g %.12g %.12g)", ea.x, ea.y, ea.z);
1162  LWDEBUGF(4, "e2 cross product == POINT(%.12g %.12g %.12g)", eb.x, eb.y, eb.z);
1163  LWDEBUGF(4, "fabs(dot_product(ea, eb)) == %.14g", fabs(dot_product(&ea, &eb)));
1164  if ( FP_EQUALS(fabs(dot_product(&ea, &eb)), 1.0) )
1165  {
1166  LWDEBUGF(4, "parallel edges found! dot_product = %.12g", dot_product(&ea, &eb));
1167  /* Parallel (maybe equal) edges! */
1168  /* Hack alert, only returning ONE end of the edge right now, most do better later. */
1169  /* Hack alert #2, returning a value of 2 to indicate a co-linear crossing event. */
1170  if ( edge_contains_point(e1, &(e2->start)) )
1171  {
1172  *g = e2->start;
1173  return 2;
1174  }
1175  if ( edge_contains_point(e1, &(e2->end)) )
1176  {
1177  *g = e2->end;
1178  return 2;
1179  }
1180  if ( edge_contains_point(e2, &(e1->start)) )
1181  {
1182  *g = e1->start;
1183  return 2;
1184  }
1185  if ( edge_contains_point(e2, &(e1->end)) )
1186  {
1187  *g = e1->end;
1188  return 2;
1189  }
1190  }
1191  unit_normal(&ea, &eb, &v);
1192  LWDEBUGF(4, "v == POINT(%.12g %.12g %.12g)", v.x, v.y, v.z);
1193  g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
1194  g->lon = atan2(v.y, v.x);
1195  LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
1196  LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
1197  if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
1198  {
1199  return LW_TRUE;
1200  }
1201  else
1202  {
1203  LWDEBUG(4, "flipping point to other side of sphere");
1204  g->lat = -1.0 * g->lat;
1205  g->lon = g->lon + M_PI;
1206  if ( g->lon > M_PI )
1207  {
1208  g->lon = -1.0 * (2.0 * M_PI - g->lon);
1209  }
1210  if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
1211  {
1212  return LW_TRUE;
1213  }
1214  }
1215  return LW_FALSE;
1216 }
#define LW_FALSE
Definition: liblwgeom.h:77
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
#define FP_EQUALS(A, B)
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:615
void robust_cross_product(const GEOGRAPHIC_POINT *p, const GEOGRAPHIC_POINT *q, POINT3D *a)
Computes the cross product of two vectors using their lat, lng representations.
Definition: lwgeodetic.c:634
void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
Calculates the unit normal to two vectors, trying to avoid problems with over-narrow or over-wide cas...
Definition: lwgeodetic.c:541
static double dot_product(const POINT3D *p1, const POINT3D *p2)
Convert cartesian coordinates on unit sphere to lon/lat coordinates static void cart2ll(const POINT3D...
Definition: lwgeodetic.c:446
int edge_contains_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
Returns true if the point p is on the minor edge defined by the end points of e.
Definition: lwgeodetic.c:1034
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:170
#define rad2deg(r)
Definition: lwgeodetic.h:80
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:63
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:64
double z
Definition: liblwgeom.h:343
double x
Definition: liblwgeom.h:343
double y
Definition: liblwgeom.h:343

References dot_product(), edge_contains_point(), GEOGRAPHIC_EDGE::end, FP_EQUALS, geographic_point_equals(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_FALSE, LW_TRUE, LWDEBUG, LWDEBUGF, normalize(), rad2deg, robust_cross_product(), GEOGRAPHIC_EDGE::start, unit_normal(), POINT3D::x, POINT3D::y, and POINT3D::z.

Referenced by circ_tree_distance_tree_internal(), and test_edge_intersection().

Here is the call graph for this function:
Here is the caller graph for this function: