PostGIS  3.7.0dev-r@@SVN_REVISION@@

◆ edge_point_in_cone()

int edge_point_in_cone ( const GEOGRAPHIC_EDGE e,
const GEOGRAPHIC_POINT p 
)

Returns true if the point p is inside the cone defined by the two ends of the edge e.

Definition at line 736 of file lwgeodetic.c.

737 {
738  POINT3D vcp, vs, ve, vp;
739  double vs_dot_vcp, vp_dot_vcp;
740  geog2cart(&(e->start), &vs);
741  geog2cart(&(e->end), &ve);
742  /* Antipodal case, everything is inside. */
743  if ( vs.x == -1.0 * ve.x && vs.y == -1.0 * ve.y && vs.z == -1.0 * ve.z )
744  return LW_TRUE;
745  geog2cart(p, &vp);
746  /* The normalized sum bisects the angle between start and end. */
747  vector_sum(&vs, &ve, &vcp);
748  normalize(&vcp);
749  /* The projection of start onto the center defines the minimum similarity */
750  vs_dot_vcp = dot_product(&vs, &vcp);
751  LWDEBUGF(4,"vs_dot_vcp %.19g",vs_dot_vcp);
752  /* The projection of candidate p onto the center */
753  vp_dot_vcp = dot_product(&vp, &vcp);
754  LWDEBUGF(4,"vp_dot_vcp %.19g",vp_dot_vcp);
755  /* If p is more similar than start then p is inside the cone */
756  LWDEBUGF(4,"fabs(vp_dot_vcp - vs_dot_vcp) %.39g",fabs(vp_dot_vcp - vs_dot_vcp));
757 
758  /*
759  ** We want to test that vp_dot_vcp is >= vs_dot_vcp but there are
760  ** numerical stability issues for values that are very very nearly
761  ** equal. Unfortunately there are also values of vp_dot_vcp that are legitimately
762  ** very close to but still less than vs_dot_vcp which we also need to catch.
763  ** The tolerance of 10-17 seems to do the trick on 32-bit and 64-bit architectures,
764  ** for the test cases here.
765  ** However, tuning the tolerance value feels like a dangerous hack.
766  ** Fundamentally, the problem is that this test is so sensitive.
767  */
768 
769  /* 1.1102230246251565404236316680908203125e-16 */
770 
771  if ( vp_dot_vcp > vs_dot_vcp || fabs(vp_dot_vcp - vs_dot_vcp) < 2e-16 )
772  {
773  LWDEBUG(4, "point is in cone");
774  return LW_TRUE;
775  }
776  LWDEBUG(4, "point is not in cone");
777  return LW_FALSE;
778 }
#define LW_FALSE
Definition: liblwgeom.h:94
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:93
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:615
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
void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the sum of two vectors.
Definition: lwgeodetic.c:465
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:404
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:101
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:106
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:64
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:65
double z
Definition: liblwgeom.h:402
double x
Definition: liblwgeom.h:402
double y
Definition: liblwgeom.h:402

References dot_product(), GEOGRAPHIC_EDGE::end, geog2cart(), LW_FALSE, LW_TRUE, LWDEBUG, LWDEBUGF, normalize(), GEOGRAPHIC_EDGE::start, vector_sum(), POINT3D::x, POINT3D::y, and POINT3D::z.

Referenced by edge_contains_point(), and edge_distance_to_point().

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