PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ circ_tree_distance_tree_internal()

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 
)

Definition at line 677 of file lwgeodetic_tree.c.

678{
679 double max;
680 double d, d_min;
681 uint32_t i;
682
683 LWDEBUGF(4, "entered, min_dist=%.8g max_dist=%.8g, type1=%d, type2=%d", *min_dist, *max_dist, n1->geom_type, n2->geom_type);
684
685 // printf("-==-\n");
686 // circ_tree_print(n1, 0);
687 // printf("--\n");
688 // circ_tree_print(n2, 0);
689
690 /* Short circuit if we've already hit the minimum */
691 if( *min_dist < threshold || *min_dist == 0.0 )
692 return *min_dist;
693
694 /* If your minimum is greater than anyone's maximum, you can't hold the winner */
695 if( circ_node_min_distance(n1, n2) > *max_dist )
696 {
697 LWDEBUGF(4, "pruning pair %p, %p", n1, n2);
698 return FLT_MAX;
699 }
700
701 /* If your maximum is a new low, we'll use that as our new global tolerance */
702 max = circ_node_max_distance(n1, n2);
703 LWDEBUGF(5, "max %.8g", max);
704 if( max < *max_dist )
705 *max_dist = max;
706
707 /* Polygon on one side, primitive type on the other. Check for point-in-polygon */
708 /* short circuit. */
709 if ( n1->geom_type == POLYGONTYPE && n2->geom_type && ! lwtype_is_collection(n2->geom_type) )
710 {
711 POINT2D pt;
712 circ_tree_get_point(n2, &pt);
713 LWDEBUGF(4, "n1 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y);
714 if ( circ_tree_contains_point(n1, &pt, &(n1->pt_outside), 0, NULL) )
715 {
716 LWDEBUG(4, "it does");
717 *min_dist = 0.0;
718 geographic_point_init(pt.x, pt.y, closest1);
719 geographic_point_init(pt.x, pt.y, closest2);
720 return *min_dist;
721 }
722 }
723 /* Polygon on one side, primitive type on the other. Check for point-in-polygon */
724 /* short circuit. */
725 if ( n2->geom_type == POLYGONTYPE && n1->geom_type && ! lwtype_is_collection(n1->geom_type) )
726 {
727 POINT2D pt;
728 circ_tree_get_point(n1, &pt);
729 LWDEBUGF(4, "n2 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y);
730 if ( circ_tree_contains_point(n2, &pt, &(n2->pt_outside), 0, NULL) )
731 {
732 LWDEBUG(4, "it does");
733 geographic_point_init(pt.x, pt.y, closest1);
734 geographic_point_init(pt.x, pt.y, closest2);
735 *min_dist = 0.0;
736 return *min_dist;
737 }
738 }
739
740 /* Both leaf nodes, do a real distance calculation */
741 if( circ_node_is_leaf(n1) && circ_node_is_leaf(n2) )
742 {
743 double d;
744 GEOGRAPHIC_POINT close1, close2;
745 LWDEBUGF(4, "testing leaf pair [%d], [%d]", n1->edge_num, n2->edge_num);
746 /* One of the nodes is a point */
747 if ( n1->p1 == n1->p2 || n2->p1 == n2->p2 )
748 {
750 GEOGRAPHIC_POINT gp1, gp2;
751
752 /* Both nodes are points! */
753 if ( n1->p1 == n1->p2 && n2->p1 == n2->p2 )
754 {
755 geographic_point_init(n1->p1->x, n1->p1->y, &gp1);
756 geographic_point_init(n2->p1->x, n2->p1->y, &gp2);
757 close1 = gp1; close2 = gp2;
758 d = sphere_distance(&gp1, &gp2);
759 }
760 /* Node 1 is a point */
761 else if ( n1->p1 == n1->p2 )
762 {
763 geographic_point_init(n1->p1->x, n1->p1->y, &gp1);
764 geographic_point_init(n2->p1->x, n2->p1->y, &(e.start));
765 geographic_point_init(n2->p2->x, n2->p2->y, &(e.end));
766 close1 = gp1;
767 d = edge_distance_to_point(&e, &gp1, &close2);
768 }
769 /* Node 2 is a point */
770 else
771 {
772 geographic_point_init(n2->p1->x, n2->p1->y, &gp2);
773 geographic_point_init(n1->p1->x, n1->p1->y, &(e.start));
774 geographic_point_init(n1->p2->x, n1->p2->y, &(e.end));
775 close2 = gp2;
776 d = edge_distance_to_point(&e, &gp2, &close1);
777 }
778 LWDEBUGF(4, " got distance %g", d);
779 }
780 /* Both nodes are edges */
781 else
782 {
783 GEOGRAPHIC_EDGE e1, e2;
785 POINT3D A1, A2, B1, B2;
786 geographic_point_init(n1->p1->x, n1->p1->y, &(e1.start));
787 geographic_point_init(n1->p2->x, n1->p2->y, &(e1.end));
788 geographic_point_init(n2->p1->x, n2->p1->y, &(e2.start));
789 geographic_point_init(n2->p2->x, n2->p2->y, &(e2.end));
790 geog2cart(&(e1.start), &A1);
791 geog2cart(&(e1.end), &A2);
792 geog2cart(&(e2.start), &B1);
793 geog2cart(&(e2.end), &B2);
794 if ( edge_intersects(&A1, &A2, &B1, &B2) )
795 {
796 d = 0.0;
797 edge_intersection(&e1, &e2, &g);
798 close1 = close2 = g;
799 }
800 else
801 {
802 d = edge_distance_to_edge(&e1, &e2, &close1, &close2);
803 }
804 LWDEBUGF(4, "edge_distance_to_edge returned %g", d);
805 }
806 if ( d < *min_dist )
807 {
808 *min_dist = d;
809 *closest1 = close1;
810 *closest2 = close2;
811 }
812 return d;
813 }
814 else
815 {
816 d_min = FLT_MAX;
817 /* Drive the recursion into the COLLECTION types first so we end up with */
818 /* pairings of primitive geometries that can be forced into the point-in-polygon */
819 /* tests above. */
820 if ( n1->geom_type && lwtype_is_collection(n1->geom_type) )
821 {
823 for ( i = 0; i < n1->num_nodes; i++ )
824 {
825 d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
826 d_min = FP_MIN(d_min, d);
827 }
828 }
829 else if ( n2->geom_type && lwtype_is_collection(n2->geom_type) )
830 {
832 for ( i = 0; i < n2->num_nodes; i++ )
833 {
834 d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
835 d_min = FP_MIN(d_min, d);
836 }
837 }
838 else if ( ! circ_node_is_leaf(n1) )
839 {
841 for ( i = 0; i < n1->num_nodes; i++ )
842 {
843 d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
844 d_min = FP_MIN(d_min, d);
845 }
846 }
847 else if ( ! circ_node_is_leaf(n2) )
848 {
850 for ( i = 0; i < n2->num_nodes; i++ )
851 {
852 d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
853 d_min = FP_MIN(d_min, d);
854 }
855 }
856 else
857 {
858 /* Never get here */
859 }
860
861 return d_min;
862 }
863}
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition lwgeom.c:1196
#define POLYGONTYPE
Definition liblwgeom.h:104
#define FP_MIN(A, B)
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition lwgeodetic.c:180
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:896
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.
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.
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition lwgeodetic.c:404
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.
int circ_tree_get_point(const CIRC_NODE *node, POINT2D *pt)
Returns a POINT2D that is a vertex of the input shape.
static double circ_node_max_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
static int circ_node_is_leaf(const CIRC_NODE *node)
Internal nodes have their point references set to NULL.
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 void circ_internal_nodes_sort(CIRC_NODE **nodes, uint32_t num_nodes, const CIRC_NODE *target_node)
static double circ_node_min_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
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 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
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:390
double x
Definition liblwgeom.h:390
uint32_t num_nodes
POINT2D * p2
struct circ_node ** nodes
POINT2D * p1
POINT2D pt_outside
uint32_t geom_type

References circ_internal_nodes_sort(), circ_node_is_leaf(), circ_node_max_distance(), circ_node_min_distance(), circ_tree_contains_point(), circ_tree_distance_tree_internal(), 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(), circ_tree_distance_tree_internal(), geography_tree_closestpoint(), and geography_tree_shortestline().

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