PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ ptarray_contains_point_sphere()

int ptarray_contains_point_sphere ( const POINTARRAY pa,
const POINT2D pt_outside,
const POINT2D pt_to_test 
)

This routine returns LW_TRUE if the stabline joining the pt_outside and pt_to_test crosses the ring an odd number of times, or if the pt_to_test is on the ring boundary itself, returning LW_FALSE otherwise.

The pt_outside must be guaranteed to be outside the ring (use the geography_pt_outside() function to derive one in postgis, or the gbox_pt_outside() function if you don't mind burning CPU cycles building a gbox first).

Definition at line 3529 of file lwgeodetic.c.

References genraster::count, edge_intersects(), getPoint2d_p(), ll2cart(), LW_FALSE, LW_TRUE, LWDEBUGF, POINTARRAY::npoints, PIR_A_TOUCH_LEFT, PIR_A_TOUCH_RIGHT, PIR_B_TOUCH_RIGHT, PIR_COLINEAR, PIR_INTERSECTS, point3d_equals(), POINT2D::x, and POINT2D::y.

Referenced by lwpoly_covers_point2d(), test_ptarray_contains_point_sphere(), and test_ptarray_contains_point_sphere_iowa().

3530 {
3531  POINT3D S1, S2; /* Stab line end points */
3532  POINT3D E1, E2; /* Edge end points (3-space) */
3533  POINT2D p; /* Edge end points (lon/lat) */
3534  uint32_t count = 0, i, inter;
3535 
3536  /* Null input, not enough points for a ring? You ain't closed! */
3537  if ( ! pa || pa->npoints < 4 )
3538  return LW_FALSE;
3539 
3540  /* Set up our stab line */
3541  ll2cart(pt_to_test, &S1);
3542  ll2cart(pt_outside, &S2);
3543 
3544  /* Initialize first point */
3545  getPoint2d_p(pa, 0, &p);
3546  ll2cart(&p, &E1);
3547 
3548  /* Walk every edge and see if the stab line hits it */
3549  for ( i = 1; i < pa->npoints; i++ )
3550  {
3551  LWDEBUGF(4, "testing edge (%d)", i);
3552  LWDEBUGF(4, " start point == POINT(%.12g %.12g)", p.x, p.y);
3553 
3554  /* Read next point. */
3555  getPoint2d_p(pa, i, &p);
3556  ll2cart(&p, &E2);
3557 
3558  /* Skip over too-short edges. */
3559  if ( point3d_equals(&E1, &E2) )
3560  {
3561  continue;
3562  }
3563 
3564  /* Our test point is on an edge end! Point is "in ring" by our definition */
3565  if ( point3d_equals(&S1, &E1) )
3566  {
3567  return LW_TRUE;
3568  }
3569 
3570  /* Calculate relationship between stab line and edge */
3571  inter = edge_intersects(&S1, &S2, &E1, &E2);
3572 
3573  /* We have some kind of interaction... */
3574  if ( inter & PIR_INTERSECTS )
3575  {
3576  /* If the stabline is touching the edge, that implies the test point */
3577  /* is on the edge, so we're done, the point is in (on) the ring. */
3578  if ( (inter & PIR_A_TOUCH_RIGHT) || (inter & PIR_A_TOUCH_LEFT) )
3579  {
3580  return LW_TRUE;
3581  }
3582 
3583  /* It's a touching interaction, disregard all the left-side ones. */
3584  /* It's a co-linear intersection, ignore those. */
3585  if ( inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR )
3586  {
3587  /* Do nothing, to avoid double counts. */
3588  LWDEBUGF(4," edge (%d) crossed, disregarding to avoid double count", i, count);
3589  }
3590  else
3591  {
3592  /* Increment crossingn count. */
3593  count++;
3594  LWDEBUGF(4," edge (%d) crossed, count == %d", i, count);
3595  }
3596  }
3597  else
3598  {
3599  LWDEBUGF(4," edge (%d) did not cross", i);
3600  }
3601 
3602  /* Increment to next edge */
3603  E1 = E2;
3604  }
3605 
3606  LWDEBUGF(4,"final count == %d", count);
3607 
3608  /* An odd number of crossings implies containment! */
3609  if ( count % 2 )
3610  {
3611  return LW_TRUE;
3612  }
3613 
3614  return LW_FALSE;
3615 }
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:86
static int point3d_equals(const POINT3D *p1, const POINT3D *p2)
Utility function for ptarray_contains_point_sphere()
Definition: lwgeodetic.c:3399
#define PIR_A_TOUCH_RIGHT
Definition: lwgeodetic.h:84
#define PIR_A_TOUCH_LEFT
Definition: lwgeodetic.h:85
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
void ll2cart(const POINT2D *g, POINT3D *p)
Convert lon/lat coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:392
int count
Definition: genraster.py:56
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.
Definition: lwgeodetic.c:3424
double y
Definition: liblwgeom.h:330
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:338
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define PIR_COLINEAR
Definition: lwgeodetic.h:83
#define PIR_INTERSECTS
Definition: lwgeodetic.h:82
uint32_t npoints
Definition: liblwgeom.h:373
Here is the call graph for this function:
Here is the caller graph for this function: