From 440c90143c862ee9e00ad0d5fa8149e7a75a66c7 Mon Sep 17 00:00:00 2001
From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Fri, 2 Feb 2018 22:23:31 +0000
Subject: [PATCH] Apply node sorting before descent to CIRC_NODE trees for
 faster calculation of geography distances. Closes #4010

git-svn-id: http://svn.osgeo.org/postgis/trunk@16371 b70326c6-7e19-0410-871a-916f4a2858ee
---
 liblwgeom/lwgeodetic_tree.c | 49 +++++++++++++++++++++++++++++++++++++
 liblwgeom/lwgeodetic_tree.h |  1 +
 2 files changed, 50 insertions(+)

diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c
index 80a838d12..24c66762d 100644
--- a/liblwgeom/lwgeodetic_tree.c
+++ b/liblwgeom/lwgeodetic_tree.c
@@ -601,6 +601,51 @@ circ_tree_distance_tree(const CIRC_NODE* n1, const CIRC_NODE* n2, const SPHEROID
 	}
 }
 
+
+/***********************************************************************
+* Internal node sorting routine to make distance calculations faster?
+*/
+
+struct sort_node {
+	CIRC_NODE *node;
+	double d;
+};
+
+static int
+circ_nodes_sort_cmp(const void *a, const void *b)
+{
+	struct sort_node *node_a = (struct sort_node *)(a);
+	struct sort_node *node_b = (struct sort_node *)(b);
+	if (node_a->d < node_b->d) return -1;
+	else if (node_a->d > node_b->d) return 1;
+	else return 0;
+}
+
+static void
+circ_internal_nodes_sort(CIRC_NODE **nodes, uint32_t num_nodes, const CIRC_NODE *target_node)
+{
+	int i;
+	struct sort_node sort_nodes[CIRC_NODE_SIZE];
+
+	/* Copy incoming nodes into sorting array and calculate */
+	/* distance to the target node */
+	for (i = 0; i < num_nodes; i++)
+	{
+		sort_nodes[i].node = nodes[i];
+		sort_nodes[i].d = sphere_distance(&(nodes[i]->center), &(target_node->center));
+	}
+
+	/* Sort the nodes and copy the result back into the input array */
+	qsort(sort_nodes, num_nodes, sizeof(struct sort_node), circ_nodes_sort_cmp);
+	for (i = 0; i < num_nodes; i++)
+	{
+		nodes[i] = sort_nodes[i].node;
+	}
+	return;
+}
+
+/***********************************************************************/
+
 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)
 {
@@ -746,6 +791,7 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl
 		/* tests above. */
 		if ( n1->geom_type && lwtype_is_collection(n1->geom_type) )
 		{
+			circ_internal_nodes_sort(n1->nodes, n1->num_nodes, n2);
 			for ( i = 0; i < n1->num_nodes; i++ )
 			{
 				d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
@@ -754,6 +800,7 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl
 		}
 		else if ( n2->geom_type && lwtype_is_collection(n2->geom_type) )
 		{
+			circ_internal_nodes_sort(n2->nodes, n2->num_nodes, n1);
 			for ( i = 0; i < n2->num_nodes; i++ )
 			{
 				d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
@@ -762,6 +809,7 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl
 		}
 		else if ( ! circ_node_is_leaf(n1) )
 		{
+			circ_internal_nodes_sort(n1->nodes, n1->num_nodes, n2);
 			for ( i = 0; i < n1->num_nodes; i++ )
 			{
 				d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
@@ -770,6 +818,7 @@ circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, doubl
 		}
 		else if ( ! circ_node_is_leaf(n2) )
 		{
+			circ_internal_nodes_sort(n2->nodes, n2->num_nodes, n1);
 			for ( i = 0; i < n2->num_nodes; i++ )
 			{
 				d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
diff --git a/liblwgeom/lwgeodetic_tree.h b/liblwgeom/lwgeodetic_tree.h
index 53ec0b16a..e9cefee0f 100644
--- a/liblwgeom/lwgeodetic_tree.h
+++ b/liblwgeom/lwgeodetic_tree.h
@@ -41,6 +41,7 @@ typedef struct circ_node
 	struct circ_node** nodes;
 	int edge_num;
 	uint32_t geom_type;
+	double d;
 	POINT2D pt_outside;
 	POINT2D* p1;
 	POINT2D* p2;
-- 
2.40.0