PostGIS  3.7.0dev-r@@SVN_REVISION@@

◆ edge_intersects()

uint32_t edge_intersects ( const POINT3D A1,
const POINT3D A2,
const POINT3D B1,
const POINT3D B2 
)

Returns non-zero if edges A and B interact.

The type of interaction is given in the return value with the bitmask elements defined above.

Definition at line 3389 of file lwgeodetic.c.

3390 {
3391  POINT3D AN, BN, VN; /* Normals to plane A and plane B */
3392  double ab_dot;
3393  int a1_side, a2_side, b1_side, b2_side;
3394  int rv = PIR_NO_INTERACT;
3395 
3396  /* Normals to the A-plane and B-plane */
3397  unit_normal(A1, A2, &AN);
3398  unit_normal(B1, B2, &BN);
3399 
3400  /* Are A-plane and B-plane basically the same? */
3401  ab_dot = dot_product(&AN, &BN);
3402 
3403  /*
3404  * https://trac.osgeo.org/postgis/ticket/5765
3405  * Failure because the colinearity check was
3406  * triggering due to an overly loose equality
3407  * check here.
3408  * if ( FP_EQUALS(fabs(ab_dot), 1.0) )
3409  */
3410  if ( 1.0 - fabs(ab_dot) <= 10e-16 )
3411  {
3412  /* Co-linear case */
3413  if ( point_in_cone(A1, A2, B1) || point_in_cone(A1, A2, B2) ||
3414  point_in_cone(B1, B2, A1) || point_in_cone(B1, B2, A2) )
3415  {
3416  rv |= PIR_INTERSECTS;
3417  rv |= PIR_COLINEAR;
3418  }
3419  return rv;
3420  }
3421 
3422  /* What side of plane-A and plane-B do the end points */
3423  /* of A and B fall? */
3424  a1_side = dot_product_side(&BN, A1);
3425  a2_side = dot_product_side(&BN, A2);
3426  b1_side = dot_product_side(&AN, B1);
3427  b2_side = dot_product_side(&AN, B2);
3428 
3429  /* Both ends of A on the same side of plane B. */
3430  if ( a1_side == a2_side && a1_side != 0 )
3431  {
3432  /* No intersection. */
3433  return PIR_NO_INTERACT;
3434  }
3435 
3436  /* Both ends of B on the same side of plane A. */
3437  if ( b1_side == b2_side && b1_side != 0 )
3438  {
3439  /* No intersection. */
3440  return PIR_NO_INTERACT;
3441  }
3442 
3443  /* A straddles B and B straddles A, so... */
3444  if ( a1_side != a2_side && (a1_side + a2_side) == 0 &&
3445  b1_side != b2_side && (b1_side + b2_side) == 0 )
3446  {
3447  /* Have to check if intersection point is inside both arcs */
3448  unit_normal(&AN, &BN, &VN);
3449  if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
3450  {
3451  return PIR_INTERSECTS;
3452  }
3453 
3454  /* Have to check if intersection point is inside both arcs */
3455  vector_scale(&VN, -1);
3456  if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
3457  {
3458  return PIR_INTERSECTS;
3459  }
3460 
3461  return PIR_NO_INTERACT;
3462  }
3463 
3464  /* The rest are all intersects variants... */
3465  rv |= PIR_INTERSECTS;
3466 
3467  /* A touches B */
3468  if ( a1_side == 0 )
3469  {
3470  /* Touches at A1, A2 is on what side? */
3471  rv |= (a2_side < 0 ? PIR_A_TOUCH_RIGHT : PIR_A_TOUCH_LEFT);
3472  }
3473  else if ( a2_side == 0 )
3474  {
3475  /* Touches at A2, A1 is on what side? */
3476  rv |= (a1_side < 0 ? PIR_A_TOUCH_RIGHT : PIR_A_TOUCH_LEFT);
3477  }
3478 
3479  /* B touches A */
3480  if ( b1_side == 0 )
3481  {
3482  /* Touches at B1, B2 is on what side? */
3483  rv |= (b2_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
3484  }
3485  else if ( b2_side == 0 )
3486  {
3487  /* Touches at B2, B1 is on what side? */
3488  rv |= (b1_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
3489  }
3490 
3491  return rv;
3492 }
static int point_in_cone(const POINT3D *A1, const POINT3D *A2, const POINT3D *P)
Utility function for checking if P is within the cone defined by A1/A2.
Definition: lwgeodetic.c:3307
void vector_scale(POINT3D *n, double scale)
Scale a vector out by a factor.
Definition: lwgeodetic.c:487
static int dot_product_side(const POINT3D *p, const POINT3D *q)
Utility function for edge_intersects(), signum with a tolerance in determining if the value is zero.
Definition: lwgeodetic.c:3374
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
#define PIR_A_TOUCH_LEFT
Definition: lwgeodetic.h:91
#define PIR_COLINEAR
Definition: lwgeodetic.h:89
#define PIR_INTERSECTS
Definition: lwgeodetic.h:88
#define PIR_A_TOUCH_RIGHT
Definition: lwgeodetic.h:90
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:92
#define PIR_B_TOUCH_LEFT
Definition: lwgeodetic.h:93
#define PIR_NO_INTERACT
Bitmask elements for edge_intersects() return value.
Definition: lwgeodetic.h:87

References dot_product(), dot_product_side(), PIR_A_TOUCH_LEFT, PIR_A_TOUCH_RIGHT, PIR_B_TOUCH_LEFT, PIR_B_TOUCH_RIGHT, PIR_COLINEAR, PIR_INTERSECTS, PIR_NO_INTERACT, point_in_cone(), unit_normal(), and vector_scale().

Referenced by circ_tree_contains_point(), circ_tree_distance_tree_internal(), lwpoly_intersects_line(), ptarray_contains_point_sphere(), ptarray_distance_spheroid(), and test_edge_intersects().

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