PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ circ_tree_distance_tree_internal()

static double circ_tree_distance_tree_internal ( const CIRC_NODE n1,
const CIRC_NODE n2,
double  threshold,
double *  min_dist,
double *  max_dist,
GEOGRAPHIC_POINT closest1,
GEOGRAPHIC_POINT closest2 
)
static

Definition at line 650 of file lwgeodetic_tree.c.

References circ_internal_nodes_sort(), circ_node_is_leaf(), circ_node_max_distance(), circ_node_min_distance(), circ_tree_contains_point(), circ_tree_get_point(), sort_node::d, edge_distance_to_edge(), edge_distance_to_point(), edge_intersection(), edge_intersects(), circ_node::edge_num, GEOGRAPHIC_EDGE::end, FP_MIN, geog2cart(), geographic_point_init(), circ_node::geom_type, LWDEBUG, LWDEBUGF, lwtype_is_collection(), circ_node::nodes, circ_node::num_nodes, circ_node::p1, circ_node::p2, POLYGONTYPE, circ_node::pt_outside, sphere_distance(), GEOGRAPHIC_EDGE::start, POINT2D::x, and POINT2D::y.

Referenced by circ_tree_distance_tree().

651 {
652  double max;
653  double d, d_min;
654  uint32_t i;
655 
656  LWDEBUGF(4, "entered, min_dist=%.8g max_dist=%.8g, type1=%d, type2=%d", *min_dist, *max_dist, n1->geom_type, n2->geom_type);
657 /*
658  circ_tree_print(n1, 0);
659  circ_tree_print(n2, 0);
660 */
661 
662  /* Short circuit if we've already hit the minimum */
663  if( *min_dist < threshold || *min_dist == 0.0 )
664  return *min_dist;
665 
666  /* If your minimum is greater than anyone's maximum, you can't hold the winner */
667  if( circ_node_min_distance(n1, n2) > *max_dist )
668  {
669  LWDEBUGF(4, "pruning pair %p, %p", n1, n2);
670  return FLT_MAX;
671  }
672 
673  /* If your maximum is a new low, we'll use that as our new global tolerance */
674  max = circ_node_max_distance(n1, n2);
675  LWDEBUGF(5, "max %.8g", max);
676  if( max < *max_dist )
677  *max_dist = max;
678 
679  /* Polygon on one side, primitive type on the other. Check for point-in-polygon */
680  /* short circuit. */
681  if ( n1->geom_type == POLYGONTYPE && n2->geom_type && ! lwtype_is_collection(n2->geom_type) )
682  {
683  POINT2D pt;
684  circ_tree_get_point(n2, &pt);
685  LWDEBUGF(4, "n1 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y);
686  if ( circ_tree_contains_point(n1, &pt, &(n1->pt_outside), NULL) )
687  {
688  LWDEBUG(4, "it does");
689  *min_dist = 0.0;
690  geographic_point_init(pt.x, pt.y, closest1);
691  geographic_point_init(pt.x, pt.y, closest2);
692  return *min_dist;
693  }
694  }
695  /* Polygon on one side, primitive type on the other. Check for point-in-polygon */
696  /* short circuit. */
697  if ( n2->geom_type == POLYGONTYPE && n1->geom_type && ! lwtype_is_collection(n1->geom_type) )
698  {
699  POINT2D pt;
700  circ_tree_get_point(n1, &pt);
701  LWDEBUGF(4, "n2 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y);
702  if ( circ_tree_contains_point(n2, &pt, &(n2->pt_outside), NULL) )
703  {
704  LWDEBUG(4, "it does");
705  geographic_point_init(pt.x, pt.y, closest1);
706  geographic_point_init(pt.x, pt.y, closest2);
707  *min_dist = 0.0;
708  return *min_dist;
709  }
710  }
711 
712  /* Both leaf nodes, do a real distance calculation */
713  if( circ_node_is_leaf(n1) && circ_node_is_leaf(n2) )
714  {
715  double d;
716  GEOGRAPHIC_POINT close1, close2;
717  LWDEBUGF(4, "testing leaf pair [%d], [%d]", n1->edge_num, n2->edge_num);
718  /* One of the nodes is a point */
719  if ( n1->p1 == n1->p2 || n2->p1 == n2->p2 )
720  {
721  GEOGRAPHIC_EDGE e;
722  GEOGRAPHIC_POINT gp1, gp2;
723 
724  /* Both nodes are points! */
725  if ( n1->p1 == n1->p2 && n2->p1 == n2->p2 )
726  {
727  geographic_point_init(n1->p1->x, n1->p1->y, &gp1);
728  geographic_point_init(n2->p1->x, n2->p1->y, &gp2);
729  close1 = gp1; close2 = gp2;
730  d = sphere_distance(&gp1, &gp2);
731  }
732  /* Node 1 is a point */
733  else if ( n1->p1 == n1->p2 )
734  {
735  geographic_point_init(n1->p1->x, n1->p1->y, &gp1);
736  geographic_point_init(n2->p1->x, n2->p1->y, &(e.start));
737  geographic_point_init(n2->p2->x, n2->p2->y, &(e.end));
738  close1 = gp1;
739  d = edge_distance_to_point(&e, &gp1, &close2);
740  }
741  /* Node 2 is a point */
742  else
743  {
744  geographic_point_init(n2->p1->x, n2->p1->y, &gp1);
745  geographic_point_init(n1->p1->x, n1->p1->y, &(e.start));
746  geographic_point_init(n1->p2->x, n1->p2->y, &(e.end));
747  close1 = gp1;
748  d = edge_distance_to_point(&e, &gp1, &close2);
749  }
750  LWDEBUGF(4, " got distance %g", d);
751  }
752  /* Both nodes are edges */
753  else
754  {
755  GEOGRAPHIC_EDGE e1, e2;
757  POINT3D A1, A2, B1, B2;
758  geographic_point_init(n1->p1->x, n1->p1->y, &(e1.start));
759  geographic_point_init(n1->p2->x, n1->p2->y, &(e1.end));
760  geographic_point_init(n2->p1->x, n2->p1->y, &(e2.start));
761  geographic_point_init(n2->p2->x, n2->p2->y, &(e2.end));
762  geog2cart(&(e1.start), &A1);
763  geog2cart(&(e1.end), &A2);
764  geog2cart(&(e2.start), &B1);
765  geog2cart(&(e2.end), &B2);
766  if ( edge_intersects(&A1, &A2, &B1, &B2) )
767  {
768  d = 0.0;
769  edge_intersection(&e1, &e2, &g);
770  close1 = close2 = g;
771  }
772  else
773  {
774  d = edge_distance_to_edge(&e1, &e2, &close1, &close2);
775  }
776  LWDEBUGF(4, "edge_distance_to_edge returned %g", d);
777  }
778  if ( d < *min_dist )
779  {
780  *min_dist = d;
781  *closest1 = close1;
782  *closest2 = close2;
783  }
784  return d;
785  }
786  else
787  {
788  d_min = FLT_MAX;
789  /* Drive the recursion into the COLLECTION types first so we end up with */
790  /* pairings of primitive geometries that can be forced into the point-in-polygon */
791  /* tests above. */
792  if ( n1->geom_type && lwtype_is_collection(n1->geom_type) )
793  {
795  for ( i = 0; i < n1->num_nodes; i++ )
796  {
797  d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
798  d_min = FP_MIN(d_min, d);
799  }
800  }
801  else if ( n2->geom_type && lwtype_is_collection(n2->geom_type) )
802  {
804  for ( i = 0; i < n2->num_nodes; i++ )
805  {
806  d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
807  d_min = FP_MIN(d_min, d);
808  }
809  }
810  else if ( ! circ_node_is_leaf(n1) )
811  {
813  for ( i = 0; i < n1->num_nodes; i++ )
814  {
815  d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
816  d_min = FP_MIN(d_min, d);
817  }
818  }
819  else if ( ! circ_node_is_leaf(n2) )
820  {
822  for ( i = 0; i < n2->num_nodes; i++ )
823  {
824  d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
825  d_min = FP_MIN(d_min, d);
826  }
827  }
828  else
829  {
830  /* Never get here */
831  }
832 
833  return d_min;
834  }
835 }
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:917
double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
Calculate the distance between two edges.
Definition: lwgeodetic.c:1238
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:56
uint32_t num_nodes
uint32_t geom_type
int circ_tree_get_point(const CIRC_NODE *node, POINT2D *pt)
Returns a POINT2D that is a vertex of the input shape.
POINT2D * p2
#define POLYGONTYPE
Definition: liblwgeom.h:86
POINT2D * p1
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
static double circ_tree_distance_tree_internal(const CIRC_NODE *n1, const CIRC_NODE *n2, double threshold, double *min_dist, double *max_dist, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
#define FP_MIN(A, B)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:47
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
static int circ_node_is_leaf(const CIRC_NODE *node)
Internal nodes have their point references set to NULL.
static double circ_node_max_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:58
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1187
int circ_tree_contains_point(const CIRC_NODE *node, const POINT2D *pt, const POINT2D *pt_outside, int *on_boundary)
Walk the tree and count intersections between the stab line and the edges.
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition: lwgeom.c:1093
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
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:59
double y
Definition: liblwgeom.h:330
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:373
static double circ_node_min_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:171
POINT2D pt_outside
static void circ_internal_nodes_sort(CIRC_NODE **nodes, uint32_t num_nodes, const CIRC_NODE *target_node)
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
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.
Definition: lwgeodetic.c:1096
struct circ_node ** nodes
Here is the call graph for this function:
Here is the caller graph for this function: