PostGIS  3.0.6dev-r@@SVN_REVISION@@

◆ circ_tree_contains_point()

int circ_tree_contains_point ( const CIRC_NODE node,
const POINT2D pt,
const POINT2D pt_outside,
int  level,
int *  on_boundary 
)

Walk the tree and count intersections between the stab line and the edges.

odd => containment, even => no containment. KNOWN PROBLEM: Grazings (think of a sharp point, just touching the stabline) will be counted for one, which will throw off the count.

Definition at line 495 of file lwgeodetic_tree.c.

496 {
497  GEOGRAPHIC_POINT closest;
498  GEOGRAPHIC_EDGE stab_edge, edge;
499  POINT3D S1, S2, E1, E2;
500  double d;
501  uint32_t i, c;
502 
503  /* Construct a stabline edge from our "inside" to our known outside point */
504  geographic_point_init(pt->x, pt->y, &(stab_edge.start));
505  geographic_point_init(pt_outside->x, pt_outside->y, &(stab_edge.end));
506  geog2cart(&(stab_edge.start), &S1);
507  geog2cart(&(stab_edge.end), &S2);
508 
509  LWDEBUGF(3, "%*s entered", level, "");
510 
511  /*
512  * If the stabline doesn't cross within the radius of a node, there's no
513  * way it can cross.
514  */
515 
516  LWDEBUGF(3, "%*s :working on node %p, edge_num %d, radius %g, center POINT(%.12g %.12g)", level, "", node, node->edge_num, node->radius, rad2deg(node->center.lon), rad2deg(node->center.lat));
517  d = edge_distance_to_point(&stab_edge, &(node->center), &closest);
518  LWDEBUGF(3, "%*s :edge_distance_to_point=%g, node_radius=%g", level, "", d, node->radius);
519  if ( FP_LTEQ(d, node->radius) )
520  {
521  LWDEBUGF(3,"%*s :entering this branch (%p)", level, "", node);
522 
523  /* Return the crossing number of this leaf */
524  if ( circ_node_is_leaf(node) )
525  {
526  int inter;
527  LWDEBUGF(3, "%*s :leaf node calculation (edge %d)", level, "", node->edge_num);
528  geographic_point_init(node->p1->x, node->p1->y, &(edge.start));
529  geographic_point_init(node->p2->x, node->p2->y, &(edge.end));
530  geog2cart(&(edge.start), &E1);
531  geog2cart(&(edge.end), &E2);
532 
533  inter = edge_intersects(&S1, &S2, &E1, &E2);
534  LWDEBUGF(3, "%*s :inter = %d", level, "", inter);
535 
536  if ( inter & PIR_INTERSECTS )
537  {
538  LWDEBUGF(3,"%*s ::got stab line edge_intersection with this edge!", level, "");
539  /* To avoid double counting crossings-at-a-vertex, */
540  /* always ignore crossings at "lower" ends of edges*/
541  GEOGRAPHIC_POINT e1, e2;
542  cart2geog(&E1,&e1); cart2geog(&E2,&e2);
543 
544  LWDEBUGF(3,"%*s LINESTRING(%.15g %.15g,%.15g %.15g)", level, "",
545  pt->x, pt->y,
546  pt_outside->x, pt_outside->y
547  );
548 
549  LWDEBUGF(3,"%*s LINESTRING(%.15g %.15g,%.15g %.15g)", level, "",
550  rad2deg(e1.lon), rad2deg(e1.lat),
551  rad2deg(e2.lon), rad2deg(e2.lat)
552  );
553 
554  if ( inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR )
555  {
556  LWDEBUGF(3,"%*s ::rejecting stab line grazing by left-side edge", level, "");
557  return 0;
558  }
559  else
560  {
561  LWDEBUGF(3,"%*s ::accepting stab line intersection", level, "");
562  return 1;
563  }
564  }
565  else
566  {
567  LWDEBUGF(3,"%*s edge does not intersect", level, "");
568  }
569  }
570  /* Or, add up the crossing numbers of all children of this node. */
571  else
572  {
573  c = 0;
574  for ( i = 0; i < node->num_nodes; i++ )
575  {
576  LWDEBUGF(3,"%*s calling circ_tree_contains_point on child %d!", level, "", i);
577  c += circ_tree_contains_point(node->nodes[i], pt, pt_outside, level + 1, on_boundary);
578  }
579  return c % 2;
580  }
581  }
582  else
583  {
584  LWDEBUGF(3,"%*s skipping this branch (%p)", level, "", node);
585  }
586  return 0;
587 }
#define FP_LTEQ(A, B)
void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
Convert cartesian coordinates on unit sphere to spherical coordinates.
Definition: lwgeodetic.c:414
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:180
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:3541
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1218
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:404
#define rad2deg(r)
Definition: lwgeodetic.h:81
#define PIR_COLINEAR
Definition: lwgeodetic.h:89
#define PIR_INTERSECTS
Definition: lwgeodetic.h:88
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:92
static int circ_node_is_leaf(const CIRC_NODE *node)
Internal nodes have their point references set to NULL.
int circ_tree_contains_point(const CIRC_NODE *node, const POINT2D *pt, const POINT2D *pt_outside, int level, int *on_boundary)
Walk the tree and count intersections between the stab line and the edges.
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:64
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:65
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:63
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:54
double y
Definition: liblwgeom.h:376
double x
Definition: liblwgeom.h:376
uint32_t num_nodes
POINT2D * p2
struct circ_node ** nodes
double radius
POINT2D * p1
GEOGRAPHIC_POINT center

References cart2geog(), circ_node::center, circ_node_is_leaf(), edge_distance_to_point(), edge_intersects(), circ_node::edge_num, GEOGRAPHIC_EDGE::end, FP_LTEQ, geog2cart(), geographic_point_init(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LWDEBUGF, circ_node::nodes, circ_node::num_nodes, circ_node::p1, circ_node::p2, PIR_B_TOUCH_RIGHT, PIR_COLINEAR, PIR_INTERSECTS, rad2deg, circ_node::radius, GEOGRAPHIC_EDGE::start, POINT2D::x, and POINT2D::y.

Referenced by circ_tree_distance_tree_internal(), CircTreePIP(), test_tree_circ_pip(), and test_tree_circ_pip2().

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