[SCM] PostGIS branch master updated. 3.4.0rc1-883-gdd88d87cf

[git at osgeo.org]

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".

The branch, master has been updated
       via  dd88d87cfc3b7b7b7f059171fb090ee4b2db8836 (commit)
      from  f66b215f0f3a04b83b6c47e004b264c2a5403269 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit dd88d87cfc3b7b7b7f059171fb090ee4b2db8836
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Wed Jan 17 13:07:29 2024 -0800

    Replace the old cached point-in-polygon code with a new
    implementation of the same concept. More memory locality,
    and a more efficient search call stack, fewer unnecessary
    allocations. Tests out to 15% faster, at least on synthetic
    data. Removes some scruft and duplicate functionality that is
    now all in liblwgeom.

diff --git a/NEWS b/NEWS
index 68b4bfbbd..fa5958806 100644
--- a/NEWS
+++ b/NEWS
@@ -48,6 +48,7 @@ To take advantage of all SFCGAL featurs, SFCGAL 1.5.0+ is needed.
           in ST_AsGeoJson(record,..). (Jan Tojnar)
   - GH-744, Don't create docbook.css for the HTML manual,
             use style.css instead (Chris Mayo)
+  - Faster implementation of point-in-poly cached index (Paul Ramsey)
 
 
 PostGIS 3.4.0
diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in
index 3094837af..68c173cc3 100644
--- a/liblwgeom/Makefile.in
+++ b/liblwgeom/Makefile.in
@@ -96,6 +96,7 @@ SA_OBJS = \
 	lwin_encoded_polyline.o \
 	lwutil.o \
 	lwhomogenize.o \
+	intervaltree.o \
 	lwalgorithm.o \
 	lwstroke.o \
 	lwlinearreferencing.o \
@@ -160,6 +161,7 @@ SA_HEADERS = \
 	lwgeom_log.h \
 	lwgeom_sfcgal.h \
 	lwinline.h \
+	intervaltree.h \
 	lwin_wkt.h \
 	lwin_wkt_parse.h \
 	lwout_twkb.h \
diff --git a/liblwgeom/cunit/cu_tree.c b/liblwgeom/cunit/cu_tree.c
index ce46c1c4f..0491ce223 100644
--- a/liblwgeom/cunit/cu_tree.c
+++ b/liblwgeom/cunit/cu_tree.c
@@ -18,6 +18,7 @@
 #include "liblwgeom_internal.h"
 #include "lwgeodetic.h"
 #include "lwgeodetic_tree.h"
+#include "intervaltree.h"
 #include "cu_tester.h"
 
 
@@ -419,6 +420,191 @@ static void test_tree_circ_distance_threshold(void)
 
 }
 
+/***** Tests for IntervalTree *****/
+
+static void test_itree_once(
+	const char *polyWkt, double x, double y,
+	IntervalTreeResult iexp)
+{
+	LWPOINT *pt = lwpoint_make2d(SRID_DEFAULT, x, y);
+	LWGEOM *poly = lwgeom_from_wkt(polyWkt, LW_PARSER_CHECK_NONE);
+	if(!poly)
+		CU_FAIL_FATAL("unable to parse WKT");
+	IntervalTree *itree = itree_from_lwgeom(poly);
+	IntervalTreeResult irslt = itree_point_in_multipolygon(itree, pt);
+	itree_free(itree);
+	lwgeom_free(poly);
+	lwpoint_free(pt);
+	CU_ASSERT_EQUAL(irslt, iexp);
+}
+
+static void test_itree_square(void)
+{
+	const char *wktPoly =
+		"MULTIPOLYGON(((-1 -1, 0 -1, 1 -1, 1 0, 1 1, -1 1, -1 -1)))";
+
+	/* inside square */
+	test_itree_once(wktPoly, 0.5, 0.5, ITREE_INSIDE);
+	/* inside square at vertex */
+	test_itree_once(wktPoly, 0.0, 0.0, ITREE_INSIDE);
+	/* left of square */
+	test_itree_once(wktPoly, -3.0, 0.0, ITREE_OUTSIDE);
+	/* right of square */
+	test_itree_once(wktPoly, 3.0, 0.0, ITREE_OUTSIDE);
+	/* above square */
+	test_itree_once(wktPoly, 0.0, 3.0, ITREE_OUTSIDE);
+	/* on horizontal boundary */
+	test_itree_once(wktPoly, 0.0, 1.0, ITREE_BOUNDARY);
+	/* on vertical boundary */
+	test_itree_once(wktPoly, -1.0, 0.0, ITREE_BOUNDARY);
+}
+
+static void test_itree_hole(void)
+{
+	const char *wktPoly =
+		"POLYGON("
+		"(-10 -10, 0 -10, 10 -10, 10 0, 10 10, 0 10, -10 10, -10 -10),"
+		"(-5 -5, -5 0, -5 5, 5 5, 5 0, 5 -5, -5 -5))";
+
+	/* inside hole */
+	test_itree_once(wktPoly, 1.0, 1.0, ITREE_OUTSIDE);
+	/* inside hole at vertical vertex */
+	test_itree_once(wktPoly, 0.0, 0.0, ITREE_OUTSIDE);
+	/* left of hole */
+	test_itree_once(wktPoly, -7.0, 0.0, ITREE_INSIDE);
+	/* right of hole */
+	test_itree_once(wktPoly, 7.0, 0.0, ITREE_INSIDE);
+	/* left of hole at top edge */
+	test_itree_once(wktPoly, -7.0, 5.0, ITREE_INSIDE);
+	/* right of hole at bottom edge */
+	test_itree_once(wktPoly, 7.0, -5.0, ITREE_INSIDE);
+	/* above hole */
+	test_itree_once(wktPoly, 0.0, 7.0, ITREE_INSIDE);
+	/* on hole boundary */
+	test_itree_once(wktPoly, 0.0, 5.0, ITREE_BOUNDARY);
+	/* on hole corner */
+	test_itree_once(wktPoly, 5.0, 5.0, ITREE_BOUNDARY);
+}
+
+static void test_itree_hole_spike(void)
+{
+	const char *wktPoly =
+		"POLYGON("
+		"(-10 -10, 6 -10, 7 -10, 7.5 2, 8 -10, 9 -10, 10 -10, 10 10, -10 10, -10 2, -10 2, -10 -10),"
+		"(-5 -5, -5 5, 5 5, 5 -5, -5 -5))";
+
+	/* inside hole */
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+	/* inside hole */
+	test_itree_once(wktPoly, -2.0, -2.0, ITREE_OUTSIDE);
+	/* left of spike */
+	test_itree_once(wktPoly, 6.0, -2.0, ITREE_INSIDE);
+	/* right of spike */
+	test_itree_once(wktPoly, 9.0, -2.0, ITREE_INSIDE);
+	/* left of tip of spike */
+	test_itree_once(wktPoly, 6.0, 2.0, ITREE_INSIDE);
+	/* right of tip of spike */
+	test_itree_once(wktPoly, 9.0, 2.0, ITREE_INSIDE);
+	/* on spike tip */
+	test_itree_once(wktPoly, 7.5, 2.0, ITREE_BOUNDARY);
+	/* left of dupe vertex */
+	test_itree_once(wktPoly, -11, 2.0, ITREE_OUTSIDE);
+	/* right of dupe vertex */
+	test_itree_once(wktPoly, 11, 2.0, ITREE_OUTSIDE);
+}
+
+static void test_itree_multipoly_empty(void)
+{
+
+	/* outside empty */
+	const char *wktPoly = "POLYGON EMPTY";
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+
+	/* outside empty */
+	wktPoly = "MULTIPOLYGON EMPTY";
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+
+	/* outside collection of empty */
+	wktPoly = "MULTIPOLYGON(EMPTY, EMPTY, EMPTY)";
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+
+	/* mixed collection of empty and not */
+	wktPoly =
+		"MULTIPOLYGON(EMPTY,"
+		"((-10 -10, 6 -10, 7 -10, 7.5 2, 8 -10, 9 -10, 10 -10, 10 10, -10 10, -10 2, -10 2, -10 -10),"
+		"(-5 -5, -5 5, 5 5, 5 -5, -5 -5)))";
+
+	/* inside hole */
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+	/* inside hole */
+	test_itree_once(wktPoly, -2.0, -2.0, ITREE_OUTSIDE);
+	/* left of spike */
+	test_itree_once(wktPoly, 6.0, -2.0, ITREE_INSIDE);
+	/* right of spike */
+	test_itree_once(wktPoly, 9.0, -2.0, ITREE_INSIDE);
+	/* left of tip of spike */
+	test_itree_once(wktPoly, 6.0, 2.0, ITREE_INSIDE);
+	/* right of tip of spike */
+	test_itree_once(wktPoly, 9.0, 2.0, ITREE_INSIDE);
+	/* on spike tip */
+	test_itree_once(wktPoly, 7.5, 2.0, ITREE_BOUNDARY);
+	/* left of dupe vertex */
+	test_itree_once(wktPoly, -11, 2.0, ITREE_OUTSIDE);
+	/* right of dupe vertex */
+	test_itree_once(wktPoly, 11, 2.0, ITREE_OUTSIDE);
+
+	/* mixed collection of empty, empty in middle */
+	wktPoly =
+		"MULTIPOLYGON("
+		"((-10 -10, 6 -10, 7 -10, 7.5 2, 8 -10, 9 -10, 10 -10, 10 10, -10 10, -10 2, -10 2, -10 -10),"
+		"(-5 -5, -5 5, 5 5, 5 -5, -5 -5)),"
+		"EMPTY, ((-1 -1, 1 -1, 1 1, -1 1, -1 -1)))";
+
+	/* inside poly in hole */
+	test_itree_once(wktPoly, 0.0, 0.0, ITREE_INSIDE);
+	/* inside hole */
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+	/* inside hole */
+	test_itree_once(wktPoly, -2.0, -2.0, ITREE_OUTSIDE);
+	/* left of spike */
+	test_itree_once(wktPoly, 6.0, -2.0, ITREE_INSIDE);
+	/* right of spike */
+	test_itree_once(wktPoly, 9.0, -2.0, ITREE_INSIDE);
+	/* left of tip of spike */
+	test_itree_once(wktPoly, 6.0, 2.0, ITREE_INSIDE);
+	/* right of tip of spike */
+	test_itree_once(wktPoly, 9.0, 2.0, ITREE_INSIDE);
+	/* on spike tip */
+	test_itree_once(wktPoly, 7.5, 2.0, ITREE_BOUNDARY);
+	/* left of dupe vertex */
+	test_itree_once(wktPoly, -11, 2.0, ITREE_OUTSIDE);
+	/* right of dupe vertex */
+	test_itree_once(wktPoly, 11, 2.0, ITREE_OUTSIDE);
+}
+
+static void test_itree_degenerate_poly(void)
+{
+	/* collapsed polygon */
+	const char *wktPoly = "POLYGON((0 0, 0 0, 0 0, 0 0))";
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+
+	/* zero area polygon */
+	wktPoly = "POLYGON((0 0, 0 1, 0 1, 0 0))";
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+	test_itree_once(wktPoly, 0, 0.5, ITREE_BOUNDARY);
+	test_itree_once(wktPoly, 0, 1, ITREE_BOUNDARY);
+
+	/* zero area polygon */
+	wktPoly = "POLYGON((0 0, 0 1, 1 1, 2 2, 1 1, 0 1, 0 0))";
+	test_itree_once(wktPoly, 0.5, 0.5, ITREE_OUTSIDE);
+	test_itree_once(wktPoly, 1, 1, ITREE_BOUNDARY);
+
+	/* non finite coordinates */
+	wktPoly = "POLYGON((0 0, 0 NaN, 0 NaN, 0 0, 0 1, 0 0))";
+	test_itree_once(wktPoly, 2.0, 2.0, ITREE_OUTSIDE);
+	test_itree_once(wktPoly, 0, 0, ITREE_BOUNDARY);
+}
+
 /*
 ** Used by test harness to register the tests in this file.
 */
@@ -426,6 +612,11 @@ void tree_suite_setup(void);
 void tree_suite_setup(void)
 {
 	CU_pSuite suite = CU_add_suite("spatial_trees", NULL, NULL);
+	PG_ADD_TEST(suite, test_itree_square);
+	PG_ADD_TEST(suite, test_itree_hole);
+	PG_ADD_TEST(suite, test_itree_hole_spike);
+	PG_ADD_TEST(suite, test_itree_multipoly_empty);
+	PG_ADD_TEST(suite, test_itree_degenerate_poly);
 	PG_ADD_TEST(suite, test_tree_circ_create);
 	PG_ADD_TEST(suite, test_tree_circ_pip);
 	PG_ADD_TEST(suite, test_tree_circ_pip2);
diff --git a/liblwgeom/intervaltree.c b/liblwgeom/intervaltree.c
new file mode 100644
index 000000000..03bd52e50
--- /dev/null
+++ b/liblwgeom/intervaltree.c
@@ -0,0 +1,530 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * PostGIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * PostGIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **********************************************************************
+ *
+ * Copyright 2023 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ **********************************************************************/
+
+#include "intervaltree.h"
+
+void
+itree_free(IntervalTree *itree)
+{
+	if (itree->nodes) lwfree(itree->nodes);
+	if (itree->ringCounts) lwfree(itree->ringCounts);
+	if (itree->indexArrays)
+	{
+		for (uint32_t i = 0; i < itree->numIndexes; i++)
+		{
+			if (itree->indexArrays[i])
+				ptarray_free(itree->indexArrays[i]);
+		}
+	}
+	if (itree->indexes) lwfree(itree->indexes);
+	if (itree->indexArrays) lwfree(itree->indexArrays);
+	lwfree(itree);
+}
+
+static uint32_t
+itree_num_nodes_pointarray(const POINTARRAY *pa)
+{
+	uint32_t num_nodes = 0;
+	uint32_t level_nodes = 0;
+
+	/* not a closed polygon */
+	if (!pa || pa->npoints < 4) return 0;
+
+	/* one leaf node per edge */
+	level_nodes = pa->npoints - 1;
+
+	while(level_nodes > 1)
+	{
+		uint32_t next_level_nodes = level_nodes / ITREE_MAX_NODES;
+		next_level_nodes += ((level_nodes % ITREE_MAX_NODES) > 0);
+		num_nodes += level_nodes;
+		level_nodes = next_level_nodes;
+	}
+	/* and the root node */
+	return num_nodes + 1;
+}
+
+static uint32_t
+itree_num_nodes_polygon(const LWPOLY *poly)
+{
+	uint32_t num_nodes = 0;
+	for (uint32_t i = 0; i < poly->nrings; i++)
+	{
+		const POINTARRAY *pa = poly->rings[i];
+		num_nodes += itree_num_nodes_pointarray(pa);
+	}
+	return num_nodes;
+}
+
+static uint32_t
+itree_num_nodes_multipolygon(const LWMPOLY *mpoly)
+{
+	uint32_t num_nodes = 0;
+	for (uint32_t i = 0; i < mpoly->ngeoms; i++)
+	{
+		const LWPOLY *poly = mpoly->geoms[i];
+		num_nodes += itree_num_nodes_polygon(poly);
+	}
+	return num_nodes;
+}
+
+static uint32_t
+itree_num_rings(const LWMPOLY *mpoly)
+{
+	return lwgeom_count_rings((const LWGEOM*)mpoly);
+}
+
+static IntervalTreeNode *
+itree_new_node(IntervalTree *itree)
+{
+	IntervalTreeNode *node = NULL;
+	if (itree->numNodes >= itree->maxNodes)
+		lwerror("%s ran out of node space", __func__);
+
+	node = &(itree->nodes[itree->numNodes++]);
+	node->min = FLT_MAX;
+	node->max = FLT_MIN;
+	node->edgeIndex = UINT32_MAX;
+	node->numChildren = 0;
+	memset(node->children, 0, sizeof(node->children));
+	return node;
+}
+
+
+static uint32_t // nodes_remaining for after latest merge, stop at 1
+itree_merge_nodes(IntervalTree *itree, uint32_t nodes_remaining)
+ {
+	/*
+	 * store the starting state of the tree which gives
+	 * us the array of nodes we are going to have to merge
+	 */
+	uint32_t end_node = itree->numNodes;
+	uint32_t start_node = end_node - nodes_remaining;
+
+	/*
+	 * one parent for every ITREE_MAX_NODES children
+	 * plus one for the remainder
+	 */
+	uint32_t num_parents = nodes_remaining / ITREE_MAX_NODES;
+	num_parents += ((nodes_remaining % ITREE_MAX_NODES) > 0);
+
+	/*
+	 * each parent is composed by merging four adjacent children
+	 * this generates a useful index structure in O(n) time
+	 * because the edges are from a ring, and thus are
+	 * spatially autocorrelated, so the pre-sorting step of
+	 * building a packed index is already done for us before we even start
+	 */
+	for (uint32_t i = 0; i < num_parents; i++)
+	{
+		uint32_t children_start = start_node + i * ITREE_MAX_NODES;
+		uint32_t children_end = start_node + (i+1) * ITREE_MAX_NODES;
+		children_end = children_end > end_node ? end_node : children_end;
+
+		/*
+		 * put pointers to the children we are merging onto
+		 * the new parent node
+		 */
+		IntervalTreeNode *parent_node = itree_new_node(itree);
+		for (uint32_t j = children_start; j < children_end; j++)
+		{
+			IntervalTreeNode *child_node = &(itree->nodes[j]);
+			parent_node->min = FP_MIN(child_node->min, parent_node->min);
+			parent_node->max = FP_MAX(child_node->max, parent_node->max);
+			parent_node->edgeIndex = FP_MIN(child_node->edgeIndex, parent_node->edgeIndex);
+			parent_node->children[parent_node->numChildren++] = child_node;
+		}
+	}
+
+	/*
+	 * keep going until num_parents gets down to one and
+	 * we are at the root of the tree
+	 */
+	return num_parents;
+}
+
+static int
+itree_edge_invalid(const POINT2D *pt1, const POINT2D *pt2)
+{
+	/* zero length */
+	if (pt1->x == pt2->x && pt1->y == pt2->y)
+		return 1;
+
+	/* nan/inf coordinates */
+	if (isfinite(pt1->x) &&
+	    isfinite(pt1->y) &&
+	    isfinite(pt2->x) &&
+	    isfinite(pt2->y))
+		return 0;
+
+	return 1;
+}
+
+static void
+itree_add_pointarray(IntervalTree *itree, const POINTARRAY *pa)
+{
+	uint32_t nodes_remaining = 0;
+	uint32_t leaf_nodes = 0;
+	IntervalTreeNode *root = NULL;
+
+	/* EMPTY/unusable ring */
+	if (!pa || pa->npoints < 4)
+		lwerror("%s called with unusable ring", __func__);
+
+	/* fill in the leaf nodes */
+	for (uint32_t i = 0; i < pa->npoints-1; i++)
+	{
+		const POINT2D *pt1 = getPoint2d_cp(pa, i);
+		const POINT2D *pt2 = getPoint2d_cp(pa, i+1);
+
+		/* Do not add nodes for zero length segments */
+		if (itree_edge_invalid(pt1, pt2))
+			continue;
+
+		/* get a fresh node for each segment of the ring */
+		IntervalTreeNode *node = itree_new_node(itree);
+		node->min = FP_MIN(pt1->y, pt2->y);
+		node->max = FP_MAX(pt1->y, pt2->y);
+		node->edgeIndex = i;
+		leaf_nodes++;
+	}
+
+	/* merge leaf nodes up to parents */
+	nodes_remaining = leaf_nodes;
+	while (nodes_remaining > 1)
+		nodes_remaining = itree_merge_nodes(itree, nodes_remaining);
+
+	/* final parent is the root */
+	if (leaf_nodes > 0)
+		root = &(itree->nodes[itree->numNodes - 1]);
+	else
+		root = NULL;
+
+	/*
+	 * take a copy of the point array we built this
+	 * tree on top of so we can reference it to get
+	 * segment information later
+	 */
+	itree->indexes[itree->numIndexes] = root;
+	itree->indexArrays[itree->numIndexes] = ptarray_clone(pa);
+	itree->numIndexes += 1;
+
+	return;
+}
+
+
+static void
+itree_add_polygon(IntervalTree *itree, const LWPOLY *poly)
+{
+	if (poly->nrings == 0) return;
+
+	itree->maxNodes = itree_num_nodes_polygon(poly);
+	itree->nodes = lwalloc0(itree->maxNodes * sizeof(IntervalTreeNode));
+	itree->numNodes = 0;
+
+	itree->ringCounts = lwalloc0(sizeof(uint32_t));
+	itree->indexes = lwalloc0(poly->nrings * sizeof(IntervalTreeNode*));
+	itree->indexArrays = lwalloc0(poly->nrings * sizeof(POINTARRAY*));
+
+	for (uint32_t j = 0; j < poly->nrings; j++)
+	{
+		const POINTARRAY *pa = poly->rings[j];
+
+		/* skip empty/unclosed/invalid rings */
+		if (!pa || pa->npoints < 4)
+			continue;
+
+		itree_add_pointarray(itree, pa);
+
+		itree->ringCounts[itree->numPolys] += 1;
+	}
+	itree->numPolys = 1;
+	return;
+}
+
+
+static void
+itree_add_multipolygon(IntervalTree *itree, const LWMPOLY *mpoly)
+{
+	if (mpoly->ngeoms == 0) return;
+
+	itree->maxNodes = itree_num_nodes_multipolygon(mpoly);
+	itree->nodes = lwalloc0(itree->maxNodes * sizeof(IntervalTreeNode));
+	itree->numNodes = 0;
+
+	itree->ringCounts = lwalloc0(mpoly->ngeoms * sizeof(uint32_t));
+	itree->indexes = lwalloc0(itree_num_rings(mpoly) * sizeof(IntervalTreeNode*));
+	itree->indexArrays = lwalloc0(itree_num_rings(mpoly) * sizeof(POINTARRAY*));
+
+	for (uint32_t i = 0; i < mpoly->ngeoms; i++)
+	{
+		const LWPOLY *poly = mpoly->geoms[i];
+
+		/* skip empty polygons */
+		if (! poly || lwpoly_is_empty(poly))
+			continue;
+
+		for (uint32_t j = 0; j < poly->nrings; j++)
+		{
+			const POINTARRAY *pa = poly->rings[j];
+
+			/* skip empty/unclosed/invalid rings */
+			if (!pa || pa->npoints < 4)
+				continue;
+
+			itree_add_pointarray(itree, pa);
+
+			itree->ringCounts[itree->numPolys] += 1;
+		}
+
+		itree->numPolys += 1;
+	}
+	return;
+}
+
+
+IntervalTree *
+itree_from_lwgeom(const LWGEOM *geom)
+{
+	if (!geom) lwerror("%s called with null geometry", __func__);
+	if (lwgeom_get_type(geom) == MULTIPOLYGONTYPE)
+	{
+		IntervalTree *itree = lwalloc0(sizeof(IntervalTree));
+		itree_add_multipolygon(itree, lwgeom_as_lwmpoly(geom));
+		return itree;
+	}
+	else if (lwgeom_get_type(geom) == POLYGONTYPE)
+	{
+		IntervalTree *itree = lwalloc0(sizeof(IntervalTree));
+		itree_add_polygon(itree, lwgeom_as_lwpoly(geom));
+		return itree;
+	}
+	else
+	{
+		lwerror("%s got asked to build index on non-polygon", __func__);
+		return NULL;
+	}
+}
+
+
+/*******************************************************************************
+ * The following is based on the "Fast Winding Number Inclusion of a Point
+ * in a Polygon" algorithm by Dan Sunday.
+ * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
+ *
+ * returns: >0 for a point to the left of the segment,
+ *          <0 for a point to the right of the segment,
+ *          0 for a point on the segment
+ */
+static inline double
+itree_segment_side(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
+{
+	return ((seg2->x - seg1->x) * (point->y - seg1->y) -
+		    (point->x - seg1->x) * (seg2->y - seg1->y));
+}
+
+/*
+ * This function doesn't test that the point falls on the line defined by
+ * the two points.  It assumes that that has already been determined
+ * by having itree_segment_side return within the tolerance.  It simply checks
+ * that if the point is on the line, it is within the endpoints.
+ *
+ * returns: 1 if the point is inside the segment bounds
+ *          0 if the point is outside the segment bounds
+ */
+static int
+itree_point_on_segment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
+{
+	double maxX = FP_MAX(seg1->x, seg2->x);
+	double maxY = FP_MAX(seg1->y, seg2->y);
+	double minX = FP_MIN(seg1->x, seg2->x);
+	double minY = FP_MIN(seg1->y, seg2->y);
+
+	return point->x >= minX && point->x <= maxX &&
+	       point->y >= minY && point->y <= maxY;
+}
+
+
+static IntervalTreeResult
+itree_point_in_ring_recursive(
+	const IntervalTreeNode *node,
+	const POINTARRAY *pa,
+	const POINT2D *pt,
+	int *winding_number)
+{
+	if (!node) return ITREE_OUTSIDE;
+	/*
+	 * If Y value is not within range of node, we can
+	 * learn nothing from this node or its children, so
+	 * we exit early.
+	 */
+	uint8_t node_contains_value = FP_CONTAINS_INCL(node->min, pt->y, node->max) ? 1 : 0;
+	if (!node_contains_value)
+		return ITREE_OK;
+
+	/* This is a leaf node, so evaluate winding number */
+	if (node->numChildren == 0)
+	{
+		const POINT2D *seg1 = getPoint2d_cp(pa, node->edgeIndex);
+		const POINT2D *seg2 = getPoint2d_cp(pa, node->edgeIndex + 1);
+		double side = itree_segment_side(seg1, seg2, pt);
+
+		/* Zero length segments are ignored. */
+		// xxxx need a unit test, what about really really short segments?
+		// if (distance2d_sqr_pt_pt(seg1, seg2) < FP_EPS*FP_EPS)
+		// 	return ITREE_OK;
+
+		/* A point on the boundary of a ring is not contained. */
+		/* WAS: if (fabs(side) < 1e-12), see ticket #852 */
+		if (side == 0.0 && itree_point_on_segment(seg1, seg2, pt) == 1)
+			return ITREE_BOUNDARY;
+
+		/*
+		 * If the point is to the left of the line, and it's rising,
+		 * then the line is to the right of the point and
+		 * circling counter-clockwise, so increment.
+		 */
+		if ((seg1->y <= pt->y) && (pt->y < seg2->y) && (side > 0))
+		{
+			*winding_number = *winding_number + 1;
+		}
+
+		/*
+		 * If the point is to the right of the line, and it's falling,
+		 * then the line is to the right of the point and circling
+		 * clockwise, so decrement.
+		 */
+		else if ((seg2->y <= pt->y) && (pt->y < seg1->y) && (side < 0))
+		{
+			*winding_number = *winding_number - 1;
+		}
+
+		return ITREE_OK;
+	}
+	/* This is an interior node, so recurse downwards */
+	else
+	{
+		for (uint32_t i = 0; i < node->numChildren; i++)
+		{
+			IntervalTreeResult rslt = itree_point_in_ring_recursive(node->children[i], pa, pt, winding_number);
+			/* Short circuit and send back boundary result */
+			if (rslt == ITREE_BOUNDARY)
+				return rslt;
+		}
+	}
+	return ITREE_OK;
+}
+
+static IntervalTreeResult
+itree_point_in_ring(const IntervalTree *itree, uint32_t ringNumber, const POINT2D *pt)
+{
+	int winding_number = 0;
+	const IntervalTreeNode *node = itree->indexes[ringNumber];
+	const POINTARRAY *pa = itree->indexArrays[ringNumber];
+	IntervalTreeResult rslt = itree_point_in_ring_recursive(node, pa, pt, &winding_number);
+
+	/* Boundary case is separate from winding number */
+	if (rslt == ITREE_BOUNDARY) return rslt;
+
+	/* Not boundary, so evaluate winding number */
+	if (winding_number == 0)
+		return ITREE_OUTSIDE;
+	else
+		return ITREE_INSIDE;
+}
+
+
+/*
+ * Test if the given point falls within the given multipolygon.
+ * Assume bbox short-circuit has already been attempted.
+ * First check if point is within any of the outer rings.
+ * If not, it is outside. If so, check if the point is
+ * within the relevant inner rings. If not, it is inside.
+ */
+IntervalTreeResult
+itree_point_in_multipolygon(const IntervalTree *itree, const LWPOINT *point)
+{
+	uint32_t i = 0;
+	const POINT2D *pt;
+	IntervalTreeResult result = ITREE_OUTSIDE;
+
+	/* Empty is not within anything */
+	if (lwpoint_is_empty(point))
+		return ITREE_OUTSIDE;
+
+	/* Non-finite point is within anything */
+	pt = getPoint2d_cp(point->point, 0);
+	if (!pt || !(isfinite(pt->x) && isfinite(pt->y)))
+		return ITREE_OUTSIDE;
+
+	/* Is the point inside any of the exterior rings of the sub-polygons? */
+	for (uint32_t p = 0; p < itree->numPolys; p++)
+	{
+		uint32_t ringCount = itree->ringCounts[p];
+
+		/* Skip empty polygons */
+		if (ringCount == 0) continue;
+
+		/* Check result against exterior ring. */
+		result = itree_point_in_ring(itree, i, pt);
+
+		/* Boundary condition is a hard stop */
+		if (result == ITREE_BOUNDARY)
+			return ITREE_BOUNDARY;
+
+		/* We are inside an exterior ring! Are we outside all the holes? */
+		if (result == ITREE_INSIDE)
+		{
+			for(uint32_t r = 1; r < itree->ringCounts[p]; r++)
+			{
+				result = itree_point_in_ring(itree, i+r, pt);
+
+				/* Boundary condition is a hard stop */
+				if (result == ITREE_BOUNDARY)
+					return ITREE_BOUNDARY;
+
+				/*
+				 * Inside a hole => Outside the polygon!
+				 * But there might be other polygons lurking
+				 * inside this hole so we have to continue
+				 * and check all the exterior rings.
+				 */
+				if (result == ITREE_INSIDE)
+					goto holes_done;
+			}
+			return ITREE_INSIDE;
+		}
+
+		holes_done:
+			/* Move to first ring of next polygon */
+			i += ringCount;
+	}
+
+	/* Not in any rings */
+	return ITREE_OUTSIDE;
+}
+
+
+
+
diff --git a/liblwgeom/intervaltree.h b/liblwgeom/intervaltree.h
new file mode 100644
index 000000000..e2b08337b
--- /dev/null
+++ b/liblwgeom/intervaltree.h
@@ -0,0 +1,76 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * PostGIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * PostGIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **********************************************************************
+ *
+ * Copyright 2023 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include "liblwgeom_internal.h"
+
+#define ITREE_MAX_NODES 4
+
+typedef enum
+{
+	ITREE_INSIDE = 1,
+	ITREE_BOUNDARY = 0, /* any boundary case stops calculations */
+	ITREE_OUTSIDE = -1,
+	ITREE_OK = 2        /* calculations may continue */
+} IntervalTreeResult;
+
+typedef struct IntervalTreeNode
+{
+	double min;
+	double max;
+	struct IntervalTreeNode *children[ITREE_MAX_NODES];
+	uint32_t edgeIndex;
+	uint32_t numChildren;
+} IntervalTreeNode;
+
+/*
+ * A multi-polygon is converted into an interval tree for
+ * searching, with one indexed ring in the index for each
+ * (non-empty, non-degenerate) ring in the input.
+ * The indexes are arranged in order they are read off
+ * the multipolygon, so (P0R0, P0R1, P1R0, P1R1) etc.
+ * The ringCounts allow hopping to the exterior rings.
+ */
+typedef struct IntervalTree
+{
+	struct IntervalTreeNode *nodes;    /* store nodes here in flat array */
+	struct IntervalTreeNode **indexes; /* store pointers to index root nodes here */
+	POINTARRAY **indexArrays;          /* store original ring data here */
+	uint32_t numIndexes;               /* utilization of the indexes/indexArrays */
+	uint32_t *ringCounts;              /* number of rings in each polygon */
+	uint32_t numPolys;                 /* number of polygons */
+	uint32_t numNodes;                 /* utilization of the nodes array */
+	uint32_t maxNodes;                 /* capacity of the nodes array */
+} IntervalTree;
+
+/*
+ * Public prototypes
+ */
+IntervalTree *itree_from_lwgeom(const LWGEOM *geom);
+void itree_free(IntervalTree *itree);
+IntervalTreeResult itree_point_in_multipolygon(const IntervalTree *itree, const LWPOINT *point);
+
+
+
diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index a0f76162f..178607155 100644
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -489,5 +489,6 @@ int gbox_contains_point2d(const GBOX *g, const POINT2D *p);
 int lwpoly_contains_point(const LWPOLY *poly, const POINT2D *pt);
 POINT4D* lwmpoint_extract_points_4d(const LWMPOINT* g, uint32_t* npoints, int* input_empty);
 char* lwstrdup(const char* a);
+void* lwalloc0(size_t sz);
 
 #endif /* _LIBLWGEOM_INTERNAL_H */
diff --git a/liblwgeom/lwutil.c b/liblwgeom/lwutil.c
index 5f63acfa8..b14bb23f2 100644
--- a/liblwgeom/lwutil.c
+++ b/liblwgeom/lwutil.c
@@ -227,14 +227,20 @@ void *
 lwalloc(size_t size)
 {
 	void *mem = lwalloc_var(size);
-	LWDEBUGF(5, "lwalloc: %d@%p", size, mem);
+	return mem;
+}
+
+void *
+lwalloc0(size_t size)
+{
+	void *mem = lwalloc_var(size);
+	memset(mem, 0, size);
 	return mem;
 }
 
 void *
 lwrealloc(void *mem, size_t size)
 {
-	LWDEBUGF(5, "lwrealloc: %d@%p", size, mem);
 	return lwrealloc_var(mem, size);
 }
 
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 2581a1f17..86454313c 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -746,8 +746,10 @@ ptarray_is_closed_z(const POINTARRAY *in)
 }
 
 /**
-* Return 1 if the point is inside the POINTARRAY, -1 if it is outside,
-* and 0 if it is on the boundary.
+* Return LW_INSIDE if the point is inside the POINTARRAY,
+* LW_OUTSIDE if it is outside, and LW_BOUNDARY if it is on
+* the boundary.
+* LW_INSIDE == 1, LW_BOUNDARY == 0, LW_OUTSIDE == -1
 */
 int
 ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
@@ -755,6 +757,12 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
 	return ptarray_contains_point_partial(pa, pt, LW_TRUE, NULL);
 }
 
+
+/*
+ * The following is based on the "Fast Winding Number Inclusion of a Point
+ * in a Polygon" algorithm by Dan Sunday.
+ * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
+ */
 int
 ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
 {
diff --git a/libpgcommon/lwgeom_cache.c b/libpgcommon/lwgeom_cache.c
index 7e013895d..95b591504 100644
--- a/libpgcommon/lwgeom_cache.c
+++ b/libpgcommon/lwgeom_cache.c
@@ -34,7 +34,7 @@
 * the following kinds of objects:
 *
 *   geometries-with-trees
-*      PreparedGeometry, RTree, CIRC_TREE, RECT_TREE
+*      PreparedGeometry, IntervalTree, CIRC_TREE, RECT_TREE
 *
 * Each GenericCache* has a type, and after that
 * some data. Similar to generic LWGEOM*. Test that
diff --git a/libpgcommon/lwgeom_cache.h b/libpgcommon/lwgeom_cache.h
index bbf3e7e18..a40512829 100644
--- a/libpgcommon/lwgeom_cache.h
+++ b/libpgcommon/lwgeom_cache.h
@@ -10,8 +10,7 @@
  *
  **********************************************************************/
 
-#ifndef LWGEOM_CACHE_H_
-#define LWGEOM_CACHE_H_ 1
+#pragma once
 
 #include "postgres.h"
 #include "fmgr.h"
@@ -20,17 +19,19 @@
 #include "lwgeodetic_tree.h"
 #include "lwgeom_pg.h"
 
-
-#define TOAST_CACHE_ENTRY 0
-#define PREP_CACHE_ENTRY 1
-#define RTREE_CACHE_ENTRY 2
-#define CIRC_CACHE_ENTRY 3
-#define RECT_CACHE_ENTRY 4
-#define SRSDESC_CACHE_ENTRY 5
-#define SRID_CACHE_ENTRY 6
+enum CacheEntryEnum {
+	TOAST_CACHE_ENTRY   = 0,
+	PREP_CACHE_ENTRY    = 1,
+	ITREE_CACHE_ENTRY   = 2,
+	CIRC_CACHE_ENTRY    = 3,
+	RECT_CACHE_ENTRY    = 4,
+	SRSDESC_CACHE_ENTRY = 5,
+	SRID_CACHE_ENTRY    = 6
+};
 
 #define NUM_CACHE_ENTRIES 7
 
+
 /* Returns the MemoryContext used to store the caches */
 MemoryContext PostgisCacheContext(FunctionCallInfo fcinfo);
 
@@ -61,7 +62,7 @@ typedef struct {
 
 /*
 * Other specific geometry cache types are the
-* RTreeGeomCache - lwgeom_rtree.h
+* IntervalTreeGeomCache - lwgeom_itree.h
 * PrepGeomCache - lwgeom_geos_prepared.h
 */
 
@@ -135,4 +136,3 @@ typedef struct {
 
 int32_t GetSRIDCacheBySRS(FunctionCallInfo fcinfo, const char *srs);
 
-#endif /* LWGEOM_CACHE_H_ */
diff --git a/postgis/Makefile.in b/postgis/Makefile.in
index 0e972da7e..c5b448b13 100644
--- a/postgis/Makefile.in
+++ b/postgis/Makefile.in
@@ -44,9 +44,6 @@ SQL_OBJS = \
 	legacy.sql \
 	legacy_minimal.sql
 
-GEOM_BACKEND_OBJ = lwgeom_geos.o
-
-BACKEND_OBJ=$(GEOM_BACKEND_OBJ)
 DATA_built=$(SQL_built)
 
 ifeq (@HAVE_SPGIST@,yes)
@@ -80,12 +77,12 @@ PG_OBJS= \
 	lwgeom_spheroid.o \
 	lwgeom_ogc.o \
 	lwgeom_functions_analytic.o \
-	lwgeom_inout.o \
 	lwgeom_functions_basic.o \
+	lwgeom_inout.o \
 	lwgeom_btree.o \
 	lwgeom_box.o \
 	lwgeom_box3d.o \
-	$(BACKEND_OBJ) \
+	lwgeom_geos.o \
 	lwgeom_geos_prepared.o \
 	lwgeom_geos_clean.o \
 	lwgeom_geos_relatematch.o \
@@ -104,9 +101,9 @@ PG_OBJS= \
 	lwgeom_functions_lrs.o \
 	lwgeom_functions_temporal.o \
 	lwgeom_rectree.o \
+	lwgeom_itree.o \
 	long_xact.o \
 	lwgeom_sqlmm.o \
-	lwgeom_rtree.o \
 	lwgeom_transform.o \
 	lwgeom_window.o \
 	gserialized_typmod.o \
diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c
index e1f94b6c9..abcc222c6 100644
--- a/postgis/lwgeom_functions_analytic.c
+++ b/postgis/lwgeom_functions_analytic.c
@@ -30,8 +30,7 @@
 #include "liblwgeom_internal.h"  /* For FP comparators. */
 #include "lwgeom_pg.h"
 #include "math.h"
-#include "lwgeom_rtree.h"
-#include "lwgeom_functions_analytic.h"
+#include "lwgeom_itree.h"
 
 #include "access/htup_details.h"
 
@@ -47,85 +46,6 @@ Datum ST_IsPolygonCCW(PG_FUNCTION_ARGS);
 Datum ST_IsPolygonCW(PG_FUNCTION_ARGS);
 
 
-/*
- * return -1 iff point is outside ring pts
- * return 1 iff point is inside ring pts
- * return 0 iff point is on ring pts
- */
-static int point_in_ring(POINTARRAY *pts, const POINT2D *point)
-{
-	int wn = 0;
-	uint32_t i;
-	double side;
-	const POINT2D* seg1;
-	const POINT2D* seg2;
-
-	POSTGIS_DEBUG(2, "point_in_ring called.");
-
-	seg2 = getPoint2d_cp(pts, 0);
-	for (i=0; i<pts->npoints-1; i++)
-	{
-		seg1 = seg2;
-		seg2 = getPoint2d_cp(pts, i+1);
-
-		side = determineSide(seg1, seg2, point);
-
-		POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1->x, seg1->y, seg2->x, seg2->y);
-		POSTGIS_DEBUGF(3, "side result: %.8f", side);
-		POSTGIS_DEBUGF(3, "counterclockwise wrap %d, clockwise wrap %d", FP_CONTAINS_BOTTOM(seg1->y, point->y, seg2->y), FP_CONTAINS_BOTTOM(seg2->y, point->y, seg1->y));
-
-		/* zero length segments are ignored. */
-		if ((seg2->x == seg1->x) && (seg2->y == seg1->y))
-		{
-			POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
-
-			continue;
-		}
-
-		/* a point on the boundary of a ring is not contained. */
-		/* WAS: if (fabs(side) < 1e-12), see #852 */
-		if (side == 0.0)
-		{
-			if (isOnSegment(seg1, seg2, point) == 1)
-			{
-				POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
-
-				return 0;
-			}
-		}
-
-		/*
-		 * If the point is to the left of the line, and it's rising,
-		 * then the line is to the right of the point and
-		 * circling counter-clockwise, so increment.
-		 */
-		if ((seg1->y <= point->y) && (point->y < seg2->y) && (side > 0))
-		{
-			POSTGIS_DEBUG(3, "incrementing winding number.");
-
-			++wn;
-		}
-		/*
-		 * If the point is to the right of the line, and it's falling,
-		 * then the line is to the right of the point and circling
-		 * clockwise, so decrement.
-		 */
-		else if ((seg2->y <= point->y) && (point->y < seg1->y) && (side < 0))
-		{
-			POSTGIS_DEBUG(3, "decrementing winding number.");
-
-			--wn;
-		}
-	}
-
-	POSTGIS_DEBUGF(3, "winding number %d", wn);
-
-	if (wn == 0)
-		return -1;
-	return 1;
-}
-
-
 /***********************************************************************
  * Simple Douglas-Peucker line simplification.
  * No checks are done to avoid introduction of self-intersections.
@@ -731,191 +651,6 @@ Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
 
 }
 
-/*******************************************************************************
- * The following is based on the "Fast Winding Number Inclusion of a Point
- * in a Polygon" algorithm by Dan Sunday.
- * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
- ******************************************************************************/
-
-/*
- * returns: >0 for a point to the left of the segment,
- *          <0 for a point to the right of the segment,
- *          0 for a point on the segment
- */
-double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
-{
-	return ((seg2->x-seg1->x)*(point->y-seg1->y)-(point->x-seg1->x)*(seg2->y-seg1->y));
-}
-
-/*
- * This function doesn't test that the point falls on the line defined by
- * the two points.  It assumes that that has already been determined
- * by having determineSide return within the tolerance.  It simply checks
- * that if the point is on the line, it is within the endpoints.
- *
- * returns: 1 if the point is not outside the bounds of the segment
- *          0 if it is
- */
-int isOnSegment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
-{
-	double maxX;
-	double maxY;
-	double minX;
-	double minY;
-
-	if (seg1->x > seg2->x)
-	{
-		maxX = seg1->x;
-		minX = seg2->x;
-	}
-	else
-	{
-		maxX = seg2->x;
-		minX = seg1->x;
-	}
-	if (seg1->y > seg2->y)
-	{
-		maxY = seg1->y;
-		minY = seg2->y;
-	}
-	else
-	{
-		maxY = seg2->y;
-		minY = seg1->y;
-	}
-
-	POSTGIS_DEBUGF(3, "maxX minX/maxY minY: %.8f %.8f/%.8f %.8f", maxX, minX, maxY, minY);
-
-	if (maxX < point->x || minX > point->x)
-	{
-		POSTGIS_DEBUGF(3, "X value %.8f falls outside the range %.8f-%.8f", point->x, minX, maxX);
-
-		return 0;
-	}
-	else if (maxY < point->y || minY > point->y)
-	{
-		POSTGIS_DEBUGF(3, "Y value %.8f falls outside the range %.8f-%.8f", point->y, minY, maxY);
-
-		return 0;
-	}
-	return 1;
-}
-
-
-
-/*
- * return -1 iff point outside polygon
- * return 0 iff point on boundary
- * return 1 iff point inside polygon
- */
-int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
-{
-	uint32_t i;
-	int result, in_ring;
-	POINT2D pt;
-
-	POSTGIS_DEBUG(2, "point_in_polygon called.");
-
-	getPoint2d_p(point->point, 0, &pt);
-	/* assume bbox short-circuit has already been attempted */
-
-	/* everything is outside of an empty polygon */
-	if ( polygon->nrings == 0 ) return -1;
-
-	in_ring = point_in_ring(polygon->rings[0], &pt);
-	if ( in_ring == -1) /* outside the exterior ring */
-	{
-		POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
-		return -1;
-	}
-	result = in_ring;
-
-	for (i=1; i<polygon->nrings; i++)
-	{
-		in_ring = point_in_ring(polygon->rings[i], &pt);
-		if (in_ring == 1) /* inside a hole => outside the polygon */
-		{
-			POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
-			return -1;
-		}
-		if (in_ring == 0) /* on the edge of a hole */
-		{
-			POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
-			return 0;
-		}
-	}
-	return result; /* -1 = outside, 0 = boundary, 1 = inside */
-}
-
-/*
- * return -1 iff point outside multipolygon
- * return 0 iff point on multipolygon boundary
- * return 1 iff point inside multipolygon
- */
-int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point)
-{
-	uint32_t i, j;
-	int result, in_ring;
-	POINT2D pt;
-
-	POSTGIS_DEBUG(2, "point_in_polygon called.");
-
-	/* empty is not within anything */
-	if (lwpoint_is_empty(point)) return -1;
-
-	getPoint2d_p(point->point, 0, &pt);
-	/* assume bbox short-circuit has already been attempted */
-
-	result = -1;
-
-	for (j = 0; j < mpolygon->ngeoms; j++ )
-	{
-
-		LWPOLY *polygon = mpolygon->geoms[j];
-
-		/* everything is outside of an empty polygon */
-		if ( polygon->nrings == 0 ) continue;
-
-		in_ring = point_in_ring(polygon->rings[0], &pt);
-		if ( in_ring == -1) /* outside the exterior ring */
-		{
-			POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
-			continue;
-		}
-		if ( in_ring == 0 )
-		{
-			return 0;
-		}
-
-		result = in_ring;
-
-		for (i=1; i<polygon->nrings; i++)
-		{
-			in_ring = point_in_ring(polygon->rings[i], &pt);
-			if (in_ring == 1) /* inside a hole => outside the polygon */
-			{
-				POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
-				result = -1;
-				break;
-			}
-			if (in_ring == 0) /* on the edge of a hole */
-			{
-				POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
-				return 0;
-			}
-		}
-		if ( result != -1)
-		{
-			return result;
-		}
-	}
-	return result;
-}
-
-
-/*******************************************************************************
- * End of "Fast Winding Number Inclusion of a Point in a Polygon" derivative.
- ******************************************************************************/
 
 /**********************************************************************
  *
diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c
index 192d395a5..16ae6dd3d 100644
--- a/postgis/lwgeom_geos.c
+++ b/postgis/lwgeom_geos.c
@@ -37,11 +37,10 @@
 #include "access/htup_details.h"
 
 /* PostGIS */
-#include "lwgeom_functions_analytic.h" /* for point_in_polygon */
 #include "lwgeom_geos.h"
 #include "liblwgeom.h"
 #include "liblwgeom_internal.h"
-#include "lwgeom_rtree.h"
+#include "lwgeom_itree.h"
 #include "lwgeom_geos_prepared.h"
 #include "lwgeom_accum.h"
 
@@ -1784,7 +1783,7 @@ Datum overlaps(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( ! gbox_overlaps_2d(&box1, &box2) )
 		{
@@ -1836,7 +1835,7 @@ Datum contains(PG_FUNCTION_ARGS)
 	if (gserialized_is_empty(geom1) || gserialized_is_empty(geom2))
 		PG_RETURN_BOOL(false);
 
-	POSTGIS_DEBUG(3, "contains called.");
+	POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
 
 	/*
 	** short-circuit 1: if geom2 bounding box is not completely inside
@@ -1855,64 +1854,12 @@ Datum contains(PG_FUNCTION_ARGS)
 	*/
 	if (is_poly(geom1) && is_point(geom2))
 	{
-		SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
-		SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
-		const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
-		const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
-		RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
-		int retval;
-
-		POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
-		if (gserialized_get_type(gpoint) == POINTTYPE)
-		{
-			LWGEOM* point = lwgeom_from_gserialized(gpoint);
-			int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
-			lwgeom_free(point);
-
-			retval = (pip_result == 1); /* completely inside */
-		}
-		else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
-		{
-			LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
-			int found_completely_inside = LW_FALSE;
-
-			retval = LW_TRUE;
-			for (uint32_t i = 0; i < mpoint->ngeoms; i++)
-			{
-				int pip_result;
-				LWPOINT* pt = mpoint->geoms[i];
-				/* We need to find at least one point that's completely inside the
-				 * polygons (pip_result == 1).  As long as we have one point that's
-				 * completely inside, we can have as many as we want on the boundary
-				 * itself. (pip_result == 0)
-				 */
-				if (lwpoint_is_empty(pt)) continue;
-				pip_result = pip_short_circuit(cache, pt, gpoly);
-				if (pip_result == 1)
-					found_completely_inside = LW_TRUE;
-
-				if (pip_result == -1) /* completely outside */
-				{
-					retval = LW_FALSE;
-					break;
-				}
-			}
-
-			retval = retval && found_completely_inside;
-			lwmpoint_free(mpoint);
-		}
-		else
-		{
-			/* Never get here */
-			elog(ERROR,"Type isn't point or multipoint!");
-			PG_RETURN_BOOL(false);
-		}
-
-		return retval > 0;
-	}
-	else
-	{
-		POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
+		const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom2);
+		LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
+		IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom1);
+		bool result = itree_pip_contains(itree, lwpt);
+		lwgeom_free(lwpt);
+		PG_RETURN_BOOL(result);
 	}
 
 	initGEOS(lwpgnotice, lwgeom_geos_error);
@@ -1982,7 +1929,7 @@ Datum containsproperly(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( ! gbox_contains_2d(&box1, &box2) )
 			PG_RETURN_BOOL(false);
@@ -2037,6 +1984,7 @@ Datum covers(PG_FUNCTION_ARGS)
 	GBOX box1, box2;
 	PrepGeomCache *prep_cache;
 
+	POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
 
 	/* A.Covers(Empty) == FALSE */
 	if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
@@ -2049,7 +1997,7 @@ Datum covers(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( ! gbox_contains_2d(&box1, &box2) )
 		{
@@ -2062,53 +2010,12 @@ Datum covers(PG_FUNCTION_ARGS)
 	 */
 	if (is_poly(geom1) && is_point(geom2))
 	{
-		SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
-		SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
-		const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
-		const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
-		RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
-		int retval;
-
-		POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
-		if (gserialized_get_type(gpoint) == POINTTYPE)
-		{
-			LWGEOM* point = lwgeom_from_gserialized(gpoint);
-			int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
-			lwgeom_free(point);
-
-			retval = (pip_result != -1); /* not outside */
-		}
-		else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
-		{
-			LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
-			uint32_t i;
-
-			retval = LW_TRUE;
-			for (i = 0; i < mpoint->ngeoms; i++)
-			{
-				LWPOINT *pt = mpoint->geoms[i];
-				if (lwpoint_is_empty(pt)) continue;
-				if (pip_short_circuit(cache, pt, gpoly) == -1)
-				{
-					retval = LW_FALSE;
-					break;
-				}
-			}
-
-			lwmpoint_free(mpoint);
-		}
-		else
-		{
-			/* Never get here */
-			elog(ERROR,"Type isn't point or multipoint!");
-			PG_RETURN_NULL();
-		}
-
-		PG_RETURN_BOOL(retval);
-	}
-	else
-	{
-		POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
+		const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom2);
+		LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
+		IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom1);
+		bool result = itree_pip_covers(itree, lwpt);
+		lwgeom_free(lwpt);
+		PG_RETURN_BOOL(result);
 	}
 
 	initGEOS(lwpgnotice, lwgeom_geos_error);
@@ -2171,6 +2078,8 @@ Datum coveredby(PG_FUNCTION_ARGS)
 
 	gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
 
+	POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
+
 	/* A.CoveredBy(Empty) == FALSE */
 	if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
 		PG_RETURN_BOOL(false);
@@ -2180,7 +2089,7 @@ Datum coveredby(PG_FUNCTION_ARGS)
 	 * geom2 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( ! gbox_contains_2d(&box2, &box1) )
 		{
@@ -2195,53 +2104,12 @@ Datum coveredby(PG_FUNCTION_ARGS)
 	 */
 	if (is_point(geom1) && is_poly(geom2))
 	{
-		SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
-		SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
-		const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
-		const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
-		RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
-		int retval;
-
-		POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
-		if (gserialized_get_type(gpoint) == POINTTYPE)
-		{
-			LWGEOM* point = lwgeom_from_gserialized(gpoint);
-			int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
-			lwgeom_free(point);
-
-			retval = (pip_result != -1); /* not outside */
-		}
-		else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
-		{
-			LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
-			uint32_t i;
-
-			retval = LW_TRUE;
-			for (i = 0; i < mpoint->ngeoms; i++)
-			{
-				LWPOINT *pt = mpoint->geoms[i];
-				if (lwpoint_is_empty(pt)) continue;
-				if (pip_short_circuit(cache, pt, gpoly) == -1)
-				{
-					retval = LW_FALSE;
-					break;
-				}
-			}
-
-			lwmpoint_free(mpoint);
-		}
-		else
-		{
-			/* Never get here */
-			elog(ERROR,"Type isn't point or multipoint!");
-			PG_RETURN_NULL();
-		}
-
-		PG_RETURN_BOOL(retval);
-	}
-	else
-	{
-		POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
+		const GSERIALIZED *gpoint = shared_gserialized_get(shared_geom1);
+		IntervalTree *itree = GetIntervalTree(fcinfo, shared_geom2);
+		LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
+		bool result = itree_pip_covers(itree, lwpt);
+		lwgeom_free(lwpt);
+		PG_RETURN_BOOL(result);
 	}
 
 	initGEOS(lwpgnotice, lwgeom_geos_error);
@@ -2291,7 +2159,7 @@ Datum crosses(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
 		{
@@ -2347,7 +2215,7 @@ Datum ST_Intersects(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
 			PG_RETURN_BOOL(false);
@@ -2355,51 +2223,19 @@ Datum ST_Intersects(PG_FUNCTION_ARGS)
 
 	/*
 	 * short-circuit 2: if the geoms are a point and a polygon,
-	 * call the point_outside_polygon function.
+	 * call the itree_pip_intersects function.
 	 */
-	if ((is_point(geom1) && is_poly(geom2)) || (is_poly(geom1) && is_point(geom2)))
+	if ((is_point(geom1) && is_poly(geom2)) ||
+	    (is_point(geom2) && is_poly(geom1)))
 	{
 		SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
 		SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
-		const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
 		const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
-		RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
-		int retval;
-
-		POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
-		if (gserialized_get_type(gpoint) == POINTTYPE)
-		{
-			LWGEOM* point = lwgeom_from_gserialized(gpoint);
-			int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
-			lwgeom_free(point);
-
-			retval = (pip_result != -1); /* not outside */
-		}
-		else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
-		{
-			LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
-			retval = LW_FALSE;
-			for (uint32_t i = 0; i < mpoint->ngeoms; i++)
-			{
-				LWPOINT *pt = mpoint->geoms[i];
-				if (lwpoint_is_empty(pt)) continue;
-				if (pip_short_circuit(cache, pt, gpoly) != -1) /* not outside */
-				{
-					retval = LW_TRUE;
-					break;
-				}
-			}
-
-			lwmpoint_free(mpoint);
-		}
-		else
-		{
-			/* Never get here */
-			elog(ERROR,"Type isn't point or multipoint!");
-			PG_RETURN_NULL();
-		}
-
-		PG_RETURN_BOOL(retval);
+		LWGEOM *lwpt = lwgeom_from_gserialized(gpoint);
+		IntervalTree *itree = GetIntervalTree(fcinfo, shared_gpoly);
+		bool result = itree_pip_intersects(itree, lwpt);
+		lwgeom_free(lwpt);
+		PG_RETURN_BOOL(result);
 	}
 
 	initGEOS(lwpgnotice, lwgeom_geos_error);
@@ -2468,7 +2304,7 @@ Datum touches(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return FALSE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-			gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
 		{
@@ -2525,7 +2361,7 @@ Datum disjoint(PG_FUNCTION_ARGS)
 	 * geom1 bounding box we can return TRUE.
 	 */
 	if ( gserialized_get_gbox_p(geom1, &box1) &&
-	        gserialized_get_gbox_p(geom2, &box2) )
+	     gserialized_get_gbox_p(geom2, &box2) )
 	{
 		if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
 		{
diff --git a/postgis/lwgeom_geos_prepared.h b/postgis/lwgeom_geos_prepared.h
index 5f3417bf2..e945d95c5 100644
--- a/postgis/lwgeom_geos_prepared.h
+++ b/postgis/lwgeom_geos_prepared.h
@@ -22,9 +22,7 @@
  *
  **********************************************************************/
 
-
-#ifndef LWGEOM_GEOS_PREPARED_H_
-#define LWGEOM_GEOS_PREPARED_H_ 1
+#pragma once
 
 #include "postgres.h"
 #include "fmgr.h"
@@ -39,20 +37,16 @@
 #include "lwgeom_cache.h"
 
 /*
-* Cache structure. We use GSERIALIZED as keys so no transformations
-* are needed before we memcmp them with other keys. We store the
-* size to avoid having to calculate the size every time.
-* The argnum gives the number of function arguments we are caching.
-* Intersects requires that both arguments be checked for cacheability,
-* while Contains only requires that the containing argument be checked.
-* Both the Geometry and the PreparedGeometry have to be cached,
-* because the PreparedGeometry contains a reference to the geometry.
-*
-* Note that the first 6 entries are part of the common GeomCache
-* structure and have to remain in order to allow the overall caching
-* system to share code (the cache checking code is common between
-* prepared geometry, circtrees, recttrees, and rtrees).
-*/
+ * Cache structure. The common components used across all
+ * caches are in the GeomCache. That contains SHARED_GSERIALIZED
+ * and the argnum that indicates which argument we are caching
+ * prepared geometry for.
+ * The argnum gives the number of function arguments we are caching.
+ * Intersects requires that both arguments be checked for cacheability,
+ * while Contains only requires that the containing argument be checked.
+ * Both the Geometry and the PreparedGeometry have to be cached,
+ * because the PreparedGeometry contains a reference to the geometry.
+ */
 typedef struct {
 	GeomCache                   gcache;
 	MemoryContext               context_statement;
@@ -63,13 +57,12 @@ typedef struct {
 
 
 /*
-** Get the current cache, given the input geometries.
-** Function will create cache if none exists, and prepare geometries in
-** cache if necessary, or pull an existing cache if possible.
-**
-** If you are only caching one argument (e.g., in contains) supply 0 as the
-** value for pg_geom2.
-*/
+ * Get the current cache, given the input geometries.
+ * Function will create cache if none exists, and prepare geometries in
+ * cache if necessary, or pull an existing cache if possible.
+ *
+ * If you are only caching one argument (e.g., in contains) supply 0 as the
+ * value for pg_geom2.
+ */
 PrepGeomCache *GetPrepGeomCache(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *pg_geom1, SHARED_GSERIALIZED *pg_geom2);
 
-#endif /* LWGEOM_GEOS_PREPARED_H_ */
diff --git a/postgis/lwgeom_itree.c b/postgis/lwgeom_itree.c
new file mode 100644
index 000000000..046f61fee
--- /dev/null
+++ b/postgis/lwgeom_itree.c
@@ -0,0 +1,292 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * PostGIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * PostGIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **********************************************************************
+ *
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ **********************************************************************/
+
+/* PostgreSQL */
+#include "postgres.h"
+#include "funcapi.h"
+#include "fmgr.h"
+
+/* Liblwgeom */
+#include "liblwgeom.h"
+#include "liblwgeom_internal.h"  /* For FP comparators. */
+#include "intervaltree.h"
+
+/* PostGIS */
+#include "lwgeom_pg.h"
+#include "lwgeom_cache.h"
+#include "lwgeom_itree.h"
+
+
+/* Prototypes */
+Datum ST_IntersectsIntervalTree(PG_FUNCTION_ARGS);
+
+
+/**********************************************************************
+* IntervalTree Caching support
+**********************************************************************/
+
+/*
+ * Specific cache types include all the generic GeomCache slots and
+ * their own slots for their index trees. .
+ */
+typedef struct {
+	GeomCache           gcache;
+	IntervalTree        *itree;
+} IntervalTreeGeomCache;
+
+/**
+* Builder, freeer and public accessor for cached IntervalTrees
+*/
+static int
+IntervalTreeBuilder(const LWGEOM *lwgeom, GeomCache *geomcache)
+{
+	IntervalTreeGeomCache *cache = (IntervalTreeGeomCache*)geomcache;
+	IntervalTree *itree = itree_from_lwgeom(lwgeom);
+
+	if (cache->itree)
+	{
+		itree_free(cache->itree);
+		cache->itree = 0;
+	}
+	if (!itree)
+		return LW_FAILURE;
+
+	cache->itree = itree;
+	return LW_SUCCESS;
+}
+
+static int
+IntervalTreeFreer(GeomCache *geomcache)
+{
+	IntervalTreeGeomCache *cache = (IntervalTreeGeomCache*)geomcache;
+	if (cache->itree)
+	{
+		itree_free(cache->itree);
+		cache->itree = 0;
+		cache->gcache.argnum = 0;
+	}
+	return LW_SUCCESS;
+}
+
+static GeomCache *
+IntervalTreeAllocator(void)
+{
+	IntervalTreeGeomCache *cache = palloc(sizeof(IntervalTreeGeomCache));
+	memset(cache, 0, sizeof(IntervalTreeGeomCache));
+	return (GeomCache*)cache;
+}
+
+static GeomCacheMethods IntervalTreeCacheMethods =
+{
+	ITREE_CACHE_ENTRY,
+	IntervalTreeBuilder,
+	IntervalTreeFreer,
+	IntervalTreeAllocator
+};
+
+/*
+ * Always return an IntervalTree, even if we don't have one
+ * cached already so that the caller can use it to do a
+ * index assisted p-i-p every time.
+ */
+IntervalTree *
+GetIntervalTree(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1)
+{
+	const GSERIALIZED *poly1;
+	LWGEOM *lwpoly1;
+
+	IntervalTree *itree = NULL;
+	IntervalTreeGeomCache *cache = (IntervalTreeGeomCache*)GetGeomCache(fcinfo, &IntervalTreeCacheMethods, g1, NULL);
+
+	/* Found a cached tree */
+	if (cache && cache->itree)
+		return cache->itree;
+
+	/* Build one on the fly */
+	poly1 = shared_gserialized_get(g1);
+	lwpoly1 = lwgeom_from_gserialized(poly1);
+	itree = itree_from_lwgeom(lwpoly1);
+	lwgeom_free(lwpoly1);
+	return itree;
+}
+
+/*
+ * A point must be fully inside (not on boundary) of
+ * a polygon to be contained. A multipoint must have
+ * at least one fully contained member and no members
+ * outside the polygon to be contained.
+ */
+bool itree_pip_contains(const IntervalTree *itree, const LWGEOM *lwpoints)
+{
+	if (lwgeom_get_type(lwpoints) == POINTTYPE)
+	{
+		return itree_point_in_multipolygon(itree, lwgeom_as_lwpoint(lwpoints)) == ITREE_INSIDE;
+	}
+	else if (lwgeom_get_type(lwpoints) == MULTIPOINTTYPE)
+	{
+		bool found_completely_inside = false;
+		LWMPOINT *mpoint = lwgeom_as_lwmpoint(lwpoints);
+		for (uint32_t i = 0; i < mpoint->ngeoms; i++)
+		{
+			IntervalTreeResult pip_result;
+			const LWPOINT* pt = mpoint->geoms[i];
+
+			if (lwpoint_is_empty(pt))
+				continue;
+			/*
+			 * We need to find at least one point that's completely inside the
+			 * polygons (pip_result == 1).  As long as we have one point that's
+			 * completely inside, we can have as many as we want on the boundary
+			 * itself. (pip_result == 0)
+			 */
+			pip_result = itree_point_in_multipolygon(itree, pt);
+
+			if (pip_result == ITREE_INSIDE)
+				found_completely_inside = true;
+
+			if (pip_result == ITREE_OUTSIDE)
+				return false;
+		}
+		return found_completely_inside;
+	}
+	else
+	{
+		elog(ERROR, "%s got a non-point input", __func__);
+		return false;
+	}
+}
+
+
+/*
+ * If any point in the point/multipoint is outside
+ * the polygon, then the polygon does not cover the point/multipoint.
+ */
+bool itree_pip_covers(const IntervalTree *itree, const LWGEOM *lwpoints)
+{
+	if (lwgeom_get_type(lwpoints) == POINTTYPE)
+	{
+		return itree_point_in_multipolygon(itree, lwgeom_as_lwpoint(lwpoints)) != ITREE_OUTSIDE;
+	}
+	else if (lwgeom_get_type(lwpoints) == MULTIPOINTTYPE)
+	{
+		LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwpoints);
+		for (uint32_t i = 0; i < mpoint->ngeoms; i++)
+		{
+			const LWPOINT *pt = mpoint->geoms[i];
+
+			if (lwpoint_is_empty(pt))
+				continue;
+
+			if (itree_point_in_multipolygon(itree, pt) == ITREE_OUTSIDE)
+				return false;
+		}
+		return true;
+	}
+	else
+	{
+		elog(ERROR, "%s got a non-point input", __func__);
+		return false;
+	}
+}
+
+/*
+ * A.intersects(B) implies if any member of the point/multipoint
+ * is not outside, then they intersect.
+ */
+bool itree_pip_intersects(const IntervalTree *itree, const LWGEOM *lwpoints)
+{
+	if (lwgeom_get_type(lwpoints) == POINTTYPE)
+	{
+		return itree_point_in_multipolygon(itree, lwgeom_as_lwpoint(lwpoints)) != ITREE_OUTSIDE;
+	}
+	else if (lwgeom_get_type(lwpoints) == MULTIPOINTTYPE)
+	{
+		LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwpoints);
+		for (uint32_t i = 0; i < mpoint->ngeoms; i++)
+		{
+			const LWPOINT *pt = mpoint->geoms[i];
+
+			if (lwpoint_is_empty(pt))
+				continue;
+
+			if (itree_point_in_multipolygon(itree, pt) != ITREE_OUTSIDE)
+				return true;
+		}
+		return false;
+	}
+	else
+	{
+		elog(ERROR, "%s got a non-point input", __func__);
+		return false;
+	}
+}
+
+
+/**********************************************************************
+* ST_IntersectsIntervalTree
+**********************************************************************/
+
+PG_FUNCTION_INFO_V1(ST_IntersectsIntervalTree);
+Datum ST_IntersectsIntervalTree(PG_FUNCTION_ARGS)
+{
+	GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
+	GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
+	LWGEOM *lwg1, *lwg2;
+	LWPOINT *lwpt;
+	IntervalTree *itree = NULL;
+	bool isPoly1, isPoly2, isPt1, isPt2;
+
+	/* Return false on empty arguments. */
+	if (gserialized_is_empty(g1) || gserialized_is_empty(g2))
+	{
+		PG_FREE_IF_COPY(g1, 0);
+		PG_FREE_IF_COPY(g2, 1);
+		PG_RETURN_BOOL(false);
+	}
+
+	lwg1 = lwgeom_from_gserialized(g1);
+	lwg2 = lwgeom_from_gserialized(g2);
+
+	isPoly1 = lwg1->type == POLYGONTYPE || lwg1->type == MULTIPOLYGONTYPE;
+	isPoly2 = lwg2->type == POLYGONTYPE || lwg2->type == MULTIPOLYGONTYPE;
+	isPt1 = lwg1->type == POINTTYPE;
+	isPt2 = lwg2->type == POINTTYPE;
+
+	/* Two points? Get outa here... */
+	if (isPoly1 && isPt2)
+	{
+		itree = itree_from_lwgeom(lwg1);
+		lwpt = lwgeom_as_lwpoint(lwg2);
+	}
+	else if (isPoly2 && isPt1)
+	{
+		itree = itree_from_lwgeom(lwg2);
+		lwpt = lwgeom_as_lwpoint(lwg1);
+	}
+
+	if (!itree)
+		elog(ERROR, "arguments to %s must a point and a polygon", __func__);
+
+	PG_RETURN_BOOL(itree_point_in_multipolygon(itree, lwpt) != ITREE_OUTSIDE);
+}
diff --git a/postgis/lwgeom_functions_analytic.h b/postgis/lwgeom_itree.h
similarity index 59%
rename from postgis/lwgeom_functions_analytic.h
rename to postgis/lwgeom_itree.h
index 79d05ec79..e46f95315 100644
--- a/postgis/lwgeom_functions_analytic.h
+++ b/postgis/lwgeom_itree.h
@@ -18,19 +18,25 @@
  *
  **********************************************************************
  *
- * Copyright 2001-2011 Refractions Research Inc.
+ * Copyright 2024 Paul Ramsey <pramsey at cleverelephant.ca>
  *
  **********************************************************************/
 
+#pragma once
 
-#include "lwgeom_rtree.h"
+#include "liblwgeom.h"
+#include "lwgeom_cache.h"
+#include "intervaltree.h"
 
-/*
-** Public prototypes for analytic functions.
+/**
+* Checks for a cache hit against the provided geometry and returns
+* a pre-built index structure (RTREE_POLY_CACHE) if one exists. Otherwise
+* builds a new one and returns that.
 */
+IntervalTree * GetIntervalTree(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1);
+
+bool itree_pip_intersects(const IntervalTree *itree, const LWGEOM *lwpoints);
+bool itree_pip_contains(const IntervalTree *itree, const LWGEOM *lwpoints);
+bool itree_pip_covers(const IntervalTree *itree, const LWGEOM *lwpoints);
 
-int point_in_polygon(LWPOLY *polygon, LWPOINT *point);
-int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *pont);
-double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point);
-int isOnSegment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point);
 
diff --git a/postgis/lwgeom_rectree.c b/postgis/lwgeom_rectree.c
index 4927750ee..2fce4407a 100644
--- a/postgis/lwgeom_rectree.c
+++ b/postgis/lwgeom_rectree.c
@@ -174,7 +174,7 @@ Datum ST_DistanceRectTreeCached(PG_FUNCTION_ARGS)
 	if (tree_cache && tree_cache->gcache.argnum)
 	{
 		RECT_NODE *n;
-		RECT_NODE *n_cached = tree_cache->index;;
+		RECT_NODE *n_cached = tree_cache->index;
 		if (tree_cache->gcache.argnum == 1)
 		{
 			LWGEOM *lwg2 = lwgeom_from_gserialized(g2);
diff --git a/postgis/lwgeom_rtree.c b/postgis/lwgeom_rtree.c
deleted file mode 100644
index 59cc77395..000000000
--- a/postgis/lwgeom_rtree.c
+++ /dev/null
@@ -1,681 +0,0 @@
-/**********************************************************************
- *
- * PostGIS - Spatial Types for PostgreSQL
- * http://postgis.net
- *
- * PostGIS is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 2 of the License, or
- * (at your option) any later version.
- *
- * PostGIS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
- *
- **********************************************************************
- *
- * Copyright (C) 2001-2005 Refractions Research Inc.
- *
- **********************************************************************/
-
-
-#include <assert.h>
-
-#include "../postgis_config.h"
-#include "lwgeom_pg.h"
-#include "liblwgeom.h"
-#include "liblwgeom_internal.h" /* For FP comparators. */
-#include "lwgeom_cache.h"
-#include "lwgeom_rtree.h"
-#include "lwgeom_functions_analytic.h"
-
-
-/* Prototypes */
-static void RTreeFree(RTREE_NODE* root);
-
-/**
-* Allocate a fresh clean RTREE_POLY_CACHE
-*/
-static RTREE_POLY_CACHE*
-RTreeCacheCreate()
-{
-	RTREE_POLY_CACHE* result;
-	result = lwalloc(sizeof(RTREE_POLY_CACHE));
-	memset(result, 0, sizeof(RTREE_POLY_CACHE));
-	return result;
-}
-
-/**
-* Recursively frees the child nodes, the interval and the line before
-* freeing the root node.
-*/
-static void
-RTreeFree(RTREE_NODE* root)
-{
-	POSTGIS_DEBUGF(2, "RTreeFree called for %p", root);
-
-	if (root->leftNode)
-		RTreeFree(root->leftNode);
-	if (root->rightNode)
-		RTreeFree(root->rightNode);
-	if (root->segment)
-	{
-		lwline_free(root->segment);
-	}
-	lwfree(root);
-}
-
-/**
-* Free the cache object and all the sub-objects properly.
-*/
-static void
-RTreeCacheClear(RTREE_POLY_CACHE* cache)
-{
-	int g, r, i;
-	POSTGIS_DEBUGF(2, "RTreeCacheClear called for %p", cache);
-	i = 0;
-	for (g = 0; g < cache->polyCount; g++)
-	{
-		for (r = 0; r < cache->ringCounts[g]; r++)
-		{
-			RTreeFree(cache->ringIndices[i]);
-			i++;
-		}
-	}
-	lwfree(cache->ringIndices);
-	lwfree(cache->ringCounts);
-	cache->ringIndices = 0;
-	cache->ringCounts = 0;
-	cache->polyCount = 0;
-}
-
-
-/**
- * Returns 1 if min < value <= max, 0 otherwise.
-*/
-static uint32
-IntervalIsContained(const RTREE_INTERVAL* interval, double value)
-{
-	return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
-}
-
-/**
-* Creates an interval with the total extents of the two given intervals.
-*/
-static void
-RTreeMergeIntervals(
-	const RTREE_INTERVAL *inter1,
-	const RTREE_INTERVAL *inter2,
-	RTREE_INTERVAL *result)
-{
-	result->max = FP_MAX(inter1->max, inter2->max);
-	result->min = FP_MIN(inter1->min, inter2->min);
-	return;
-}
-
-
-/**
-* Creates an interior node given the children.
-*/
-static RTREE_NODE*
-RTreeCreateInteriorNode(RTREE_NODE* left, RTREE_NODE* right)
-{
-	RTREE_NODE *parent;
-
-	POSTGIS_DEBUGF(2, "RTreeCreateInteriorNode called for children %p, %p", left, right);
-
-	parent = lwalloc(sizeof(RTREE_NODE));
-	parent->leftNode = left;
-	parent->rightNode = right;
-	parent->segment = NULL;
-
-	RTreeMergeIntervals(
-		&(left->interval), &(right->interval),
-		&(parent->interval));
-
-	POSTGIS_DEBUGF(3, "RTreeCreateInteriorNode returning %p", parent);
-
-	return parent;
-}
-
-/**
-* Creates a leaf node given the pointer to the start point of the segment.
-*/
-static RTREE_NODE*
-RTreeCreateLeafNode(POINTARRAY* pa, uint32_t startPoint)
-{
-	RTREE_NODE *parent;
-	LWLINE *line;
-	double value1;
-	double value2;
-	POINT4D tmp;
-	POINTARRAY *npa;
-
-	POSTGIS_DEBUGF(2, "RTreeCreateLeafNode called for point %d of %p", startPoint, pa);
-
-	if (pa->npoints < startPoint + 2)
-	{
-		lwpgerror("RTreeCreateLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);
-	}
-
-	/*
-	 * The given point array will be part of a geometry that will be freed
-	 * independently of the index.	Since we may want to cache the index,
-	 * we must create independent arrays.
-	 */
-	npa = ptarray_construct_empty(0,0,2);
-
-	getPoint4d_p(pa, startPoint, &tmp);
-	value1 = tmp.y;
-	ptarray_append_point(npa,&tmp,LW_TRUE);
-
-	getPoint4d_p(pa, startPoint+1, &tmp);
-	value2 = tmp.y;
-	ptarray_append_point(npa,&tmp,LW_TRUE);
-
-	line = lwline_construct(SRID_UNKNOWN, NULL, npa);
-
-	parent = lwalloc(sizeof(RTREE_NODE));
-	parent->interval.min = FP_MIN(value1, value2);
-	parent->interval.max = FP_MAX(value1, value2);
-	parent->segment = line;
-	parent->leftNode = NULL;
-	parent->rightNode = NULL;
-
-	POSTGIS_DEBUGF(3, "RTreeCreateLeafNode returning %p", parent);
-
-	return parent;
-}
-
-/**
-* Creates an rtree given a pointer to the point array.
-* Must copy the point array.
-*/
-static RTREE_NODE*
-RTreeCreate(POINTARRAY* pointArray)
-{
-	RTREE_NODE* root;
-	RTREE_NODE** nodes = lwalloc(pointArray->npoints * sizeof(RTREE_NODE*));
-	uint32_t i, nodeCount;
-	uint32_t childNodes, parentNodes;
-
-	POSTGIS_DEBUGF(2, "RTreeCreate called with pointarray %p", pointArray);
-
-	nodeCount = pointArray->npoints - 1;
-
-	POSTGIS_DEBUGF(3, "Total leaf nodes: %d", nodeCount);
-
-	/*
-	 * Create a leaf node for every line segment.
-	 */
-	for (i = 0; i < nodeCount; i++)
-	{
-		nodes[i] = RTreeCreateLeafNode(pointArray, i);
-	}
-
-	/*
-	 * Next we group nodes by pairs.  If there's an odd number of nodes,
-	 * we bring the last node up a level as is.	 Continue until we have
-	 * a single top node.
-	 */
-	childNodes = nodeCount;
-	parentNodes = nodeCount / 2;
-	while (parentNodes > 0)
-	{
-		POSTGIS_DEBUGF(3, "Merging %d children into %d parents.", childNodes, parentNodes);
-
-		i = 0;
-		while (i < parentNodes)
-		{
-			nodes[i] = RTreeCreateInteriorNode(nodes[i*2], nodes[i*2+1]);
-			i++;
-		}
-		/*
-		 * Check for an odd numbered final node.
-		 */
-		if (parentNodes * 2 < childNodes)
-		{
-			POSTGIS_DEBUGF(3, "Shuffling child %d to parent %d", childNodes - 1, i);
-
-			nodes[i] = nodes[childNodes - 1];
-			parentNodes++;
-		}
-		childNodes = parentNodes;
-		parentNodes = parentNodes / 2;
-	}
-
-	root = nodes[0];
-	lwfree(nodes);
-	POSTGIS_DEBUGF(3, "RTreeCreate returning %p", root);
-
-	return root;
-}
-
-
-/**
-* Merges two multilinestrings into a single multilinestring.
-*/
-static LWMLINE*
-RTreeMergeMultiLines(LWMLINE *line1, LWMLINE *line2)
-{
-	LWGEOM **geoms;
-	LWCOLLECTION *col;
-	uint32_t i, j, ngeoms;
-
-	POSTGIS_DEBUGF(2, "RTreeMergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, line1->type, line2, line2->ngeoms, line2->type);
-
-	ngeoms = line1->ngeoms + line2->ngeoms;
-	geoms = lwalloc(sizeof(LWGEOM *) * ngeoms);
-
-	j = 0;
-	for (i = 0; i < line1->ngeoms; i++, j++)
-	{
-		geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]);
-	}
-	for (i = 0; i < line2->ngeoms; i++, j++)
-	{
-		geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]);
-	}
-	col = lwcollection_construct(MULTILINETYPE, SRID_UNKNOWN, NULL, ngeoms, geoms);
-
-	POSTGIS_DEBUGF(3, "RTreeMergeMultiLines returning %p, %d, %d", col, col->ngeoms, col->type);
-
-	return (LWMLINE *)col;
-}
-
-
-/**
-* Callback function sent into the GetGeomCache generic caching system. Given a
-* LWGEOM* this function builds and stores an RTREE_POLY_CACHE into the provided
-* GeomCache object.
-*/
-static int
-RTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
-{
-	uint32_t i, p, r;
-	LWMPOLY *mpoly;
-	LWPOLY *poly;
-	int nrings;
-	RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
-	RTREE_POLY_CACHE* currentCache;
-
-	if ( ! cache )
-		return LW_FAILURE;
-
-	if ( rtree_cache->index )
-	{
-		lwpgerror("RTreeBuilder asked to build index where one already exists.");
-		return LW_FAILURE;
-	}
-
-	if (lwgeom->type == MULTIPOLYGONTYPE)
-	{
-		POSTGIS_DEBUG(2, "RTreeBuilder MULTIPOLYGON");
-		mpoly = (LWMPOLY *)lwgeom;
-		nrings = 0;
-		/*
-		** Count the total number of rings.
-		*/
-		currentCache = RTreeCacheCreate();
-		currentCache->polyCount = mpoly->ngeoms;
-		currentCache->ringCounts = lwalloc(sizeof(int) * mpoly->ngeoms);
-		for ( i = 0; i < mpoly->ngeoms; i++ )
-		{
-			currentCache->ringCounts[i] = mpoly->geoms[i]->nrings;
-			nrings += mpoly->geoms[i]->nrings;
-		}
-		currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * nrings);
-		/*
-		** Load the array in geometry order, each outer ring followed by the inner rings
-		** associated with that outer ring
-		*/
-		i = 0;
-		for ( p = 0; p < mpoly->ngeoms; p++ )
-		{
-			for ( r = 0; r < mpoly->geoms[p]->nrings; r++ )
-			{
-				currentCache->ringIndices[i] = RTreeCreate(mpoly->geoms[p]->rings[r]);
-				i++;
-			}
-		}
-		rtree_cache->index = currentCache;
-	}
-	else if ( lwgeom->type == POLYGONTYPE )
-	{
-		POSTGIS_DEBUG(2, "RTreeBuilder POLYGON");
-		poly = (LWPOLY *)lwgeom;
-		currentCache = RTreeCacheCreate();
-		currentCache->polyCount = 1;
-		currentCache->ringCounts = lwalloc(sizeof(int));
-		currentCache->ringCounts[0] = poly->nrings;
-		/*
-		** Just load the rings on in order
-		*/
-		currentCache->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);
-		for ( i = 0; i < poly->nrings; i++ )
-		{
-			currentCache->ringIndices[i] = RTreeCreate(poly->rings[i]);
-		}
-		rtree_cache->index = currentCache;
-	}
-	else
-	{
-		/* Uh oh, shouldn't be here. */
-		lwpgerror("RTreeBuilder got asked to build index on non-polygon");
-		return LW_FAILURE;
-	}
-	return LW_SUCCESS;
-}
-
-/**
-* Callback function sent into the GetGeomCache generic caching system. On a
-* cache miss, this function clears the cached index object.
-*/
-static int
-RTreeFreer(GeomCache* cache)
-{
-	RTreeGeomCache* rtree_cache = (RTreeGeomCache*)cache;
-
-	if ( ! cache )
-		return LW_FAILURE;
-
-	if ( rtree_cache->index )
-	{
-		RTreeCacheClear(rtree_cache->index);
-		lwfree(rtree_cache->index);
-		rtree_cache->index = 0;
-		rtree_cache->gcache.argnum = 0;
-	}
-	return LW_SUCCESS;
-}
-
-static GeomCache*
-RTreeAllocator(void)
-{
-	RTreeGeomCache* cache = palloc(sizeof(RTreeGeomCache));
-	memset(cache, 0, sizeof(RTreeGeomCache));
-	return (GeomCache*)cache;
-}
-
-static GeomCacheMethods RTreeCacheMethods =
-{
-	RTREE_CACHE_ENTRY,
-	RTreeBuilder,
-	RTreeFreer,
-	RTreeAllocator
-};
-
-RTREE_POLY_CACHE *
-GetRtreeCache(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1)
-{
-	RTreeGeomCache* cache = (RTreeGeomCache*)GetGeomCache(fcinfo, &RTreeCacheMethods, g1, NULL);
-	RTREE_POLY_CACHE* index = NULL;
-
-	if ( cache )
-		index = cache->index;
-
-	return index;
-}
-
-
-/**
-* Retrieves a collection of line segments given the root and crossing value.
-* The collection is a multilinestring consisting of two point lines
-* representing the segments of the ring that may be crossed by the
-* horizontal projection line at the given y value.
-*/
-LWMLINE *RTreeFindLineSegments(RTREE_NODE *root, double value)
-{
-	LWMLINE *tmp, *result;
-	LWGEOM **lwgeoms;
-
-	POSTGIS_DEBUGF(2, "RTreeFindLineSegments called for tree %p and value %8.3f", root, value);
-
-	result = NULL;
-
-	if (!IntervalIsContained(&(root->interval), value))
-	{
-		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: not contained.", root);
-
-		return NULL;
-	}
-
-	/* If there is a segment defined for this node, include it. */
-	if (root->segment)
-	{
-		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: adding segment %p %d.", root, root->segment, root->segment->type);
-
-		lwgeoms = lwalloc(sizeof(LWGEOM *));
-		lwgeoms[0] = (LWGEOM *)root->segment;
-
-		POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", root->segment, root->segment->type, lwgeom_ndims((LWGEOM *)(root->segment)));
-
-		result = (LWMLINE *)lwcollection_construct(MULTILINETYPE, SRID_UNKNOWN, NULL, 1, lwgeoms);
-	}
-
-	/* If there is a left child node, recursively include its results. */
-	if (root->leftNode)
-	{
-		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing left.", root);
-
-		tmp = RTreeFindLineSegments(root->leftNode, value);
-		if (tmp)
-		{
-			POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, lwgeom_ndims((LWGEOM *)tmp));
-
-			if (result)
-				result = RTreeMergeMultiLines(result, tmp);
-			else
-				result = tmp;
-		}
-	}
-
-	/* Same for any right child. */
-	if (root->rightNode)
-	{
-		POSTGIS_DEBUGF(3, "RTreeFindLineSegments %p: recursing right.", root);
-
-		tmp = RTreeFindLineSegments(root->rightNode, value);
-		if (tmp)
-		{
-			POSTGIS_DEBUGF(3, "Found geom %p, type %d, dim %d", tmp, tmp->type, lwgeom_ndims((LWGEOM *)tmp));
-
-			if (result)
-				result = RTreeMergeMultiLines(result, tmp);
-			else
-				result = tmp;
-		}
-	}
-
-	return result;
-}
-
-
-
-/*
- * return -1 iff point is outside ring pts
- * return 1 iff point is inside ring pts
- * return 0 iff point is on ring pts
- */
-static int point_in_ring_rtree(RTREE_NODE *root, const POINT2D *point)
-{
-	int wn = 0;
-	uint32_t i;
-	double side;
-	const POINT2D *seg1;
-	const POINT2D *seg2;
-	LWMLINE *lines;
-
-	POSTGIS_DEBUG(2, "point_in_ring called.");
-
-	lines = RTreeFindLineSegments(root, point->y);
-	if (!lines)
-		return -1;
-
-	for (i=0; i<lines->ngeoms; i++)
-	{
-		seg1 = getPoint2d_cp(lines->geoms[i]->points, 0);
-		seg2 = getPoint2d_cp(lines->geoms[i]->points, 1);
-
-		side = determineSide(seg1, seg2, point);
-
-		POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1->x, seg1->y, seg2->x, seg2->y);
-		POSTGIS_DEBUGF(3, "side result: %.8f", side);
-		POSTGIS_DEBUGF(3, "counterclockwise wrap %d, clockwise wrap %d", FP_CONTAINS_BOTTOM(seg1->y, point->y, seg2->y), FP_CONTAINS_BOTTOM(seg2->y, point->y, seg1->y));
-
-		/* zero length segments are ignored. */
-		if (((seg2->x - seg1->x)*(seg2->x - seg1->x) + (seg2->y - seg1->y)*(seg2->y - seg1->y)) < 1e-12*1e-12)
-		{
-			POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
-			continue;
-		}
-
-		/* a point on the boundary of a ring is not contained. */
-		/* WAS: if (fabs(side) < 1e-12), see #852 */
-		if (side == 0.0)
-		{
-			if (isOnSegment(seg1, seg2, point) == 1)
-			{
-				POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
-				return 0;
-			}
-		}
-
-		/*
-		 * If the point is to the left of the line, and it's rising,
-		 * then the line is to the right of the point and
-		 * circling counter-clockwise, so increment.
-		 */
-		if ((seg1->y <= point->y) && (point->y < seg2->y) && (side > 0))
-		{
-			POSTGIS_DEBUG(3, "incrementing winding number.");
-			++wn;
-		}
-		/*
-		 * If the point is to the right of the line, and it's falling,
-		 * then the line is to the right of the point and circling
-		 * clockwise, so decrement.
-		 */
-		else if ((seg2->y <= point->y) && (point->y < seg1->y) && (side < 0))
-		{
-			POSTGIS_DEBUG(3, "decrementing winding number.");
-			--wn;
-		}
-	}
-
-	POSTGIS_DEBUGF(3, "winding number %d", wn);
-
-	if (wn == 0)
-		return -1;
-	return 1;
-}
-
-/*
- * return -1 if point outside polygon
- * return 0 if point on boundary
- * return 1 if point inside polygon
- *
- * Expected **root order is each exterior ring followed by its holes, eg. EIIEIIEI
- */
-static int
-point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point)
-{
-	int i, p, r, in_ring;
-	const POINT2D *pt;
-	int result = -1;
-
-	POSTGIS_DEBUGF(2, "point_in_multipolygon_rtree called for %p %d %p.", root, polyCount, point);
-
-	/* empty is not within anything */
-	if (lwpoint_is_empty(point)) return -1;
-
-	pt = getPoint2d_cp(point->point, 0);
-
-	/* assume bbox short-circuit has already been attempted */
-
-	i = 0; /* the current index into the root array */
-
-	/* is the point inside any of the sub-polygons? */
-	for ( p = 0; p < polyCount; p++ )
-	{
-		/* Skip empty polygons */
-		if( ringCounts[p] == 0 ) continue;
-
-		in_ring = point_in_ring_rtree(root[i], pt);
-		POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", p, in_ring);
-		if ( in_ring == -1 ) /* outside the exterior ring */
-		{
-			POSTGIS_DEBUG(3, "point_in_multipolygon_rtree: outside exterior ring.");
-		}
-		else if ( in_ring == 0 ) /* on the boundary */
-		{
-			POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of exterior ring %d", p);
-			return 0;
-		}
-		else
-		{
-			result = in_ring;
-
-			for(r = 1; r < ringCounts[p]; r++)
-			{
-				in_ring = point_in_ring_rtree(root[i+r], pt);
-				POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: interior ring (%d), point_in_ring returned %d", r, in_ring);
-				if (in_ring == 1) /* inside a hole => outside the polygon */
-				{
-					POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: within hole %d of exterior ring %d", r, p);
-					result = -1;
-					break;
-				}
-				if (in_ring == 0) /* on the edge of a hole */
-				{
-					POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of hole %d of exterior ring %d", r, p);
-					return 0;
-				}
-			}
-			/* if we have a positive result, we can short-circuit and return it */
-			if (result != -1)
-			{
-				return result;
-			}
-		}
-		/* increment the index by the total number of rings in the sub-poly */
-		/* we do this here in case we short-cutted out of the poly before looking at all the rings */
-		i += ringCounts[p];
-	}
-
-	return result; /* -1 = outside, 0 = boundary, 1 = inside */
-}
-
-
-/* Utility function that checks a LWPOINT and a GSERIALIZED poly against a cache.
- * Serialized poly may be a multipart.
- */
-int
-pip_short_circuit(RTREE_POLY_CACHE *poly_cache, LWPOINT *point, const GSERIALIZED *gpoly)
-{
-	int result;
-
-	if (poly_cache && poly_cache->ringIndices)
-	{
-		result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
-	}
-	else
-	{
-		LWGEOM* poly = lwgeom_from_gserialized(gpoly);
-		if (lwgeom_get_type(poly) == POLYGONTYPE)
-		{
-			result = point_in_polygon(lwgeom_as_lwpoly(poly), point);
-		}
-		else
-		{
-			result = point_in_multipolygon(lwgeom_as_lwmpoly(poly), point);
-		}
-		lwgeom_free(poly);
-	}
-	return result;
-}
-
diff --git a/postgis/lwgeom_rtree.h b/postgis/lwgeom_rtree.h
deleted file mode 100644
index 41db74aae..000000000
--- a/postgis/lwgeom_rtree.h
+++ /dev/null
@@ -1,87 +0,0 @@
-/**********************************************************************
- *
- * PostGIS - Spatial Types for PostgreSQL
- * http://postgis.net
- *
- * PostGIS is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 2 of the License, or
- * (at your option) any later version.
- *
- * PostGIS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
- *
- **********************************************************************
- *
- * ^copyright^
- *
- **********************************************************************/
-
-#ifndef _LWGEOM_RTREE_H
-#define _LWGEOM_RTREE_H 1
-
-#include "liblwgeom.h"
-#include "lwgeom_cache.h"
-
-/**
-* Representation for the y-axis interval spanned by an edge.
-*/
-typedef struct
-{
-	double min;
-	double max;
-}
-RTREE_INTERVAL;
-
-/**
-* The following struct and methods are used for a 1D RTree implementation,
-* described at:
-*  http://lin-ear-th-inking.blogspot.com/2007/06/packed-1-dimensional-r-tree.html
-*/
-typedef struct rtree_node
-{
-	RTREE_INTERVAL interval;
-	struct rtree_node *leftNode;
-	struct rtree_node *rightNode;
-	LWLINE *segment;
-}
-RTREE_NODE;
-
-/**
-* The tree structure used for fast P-i-P tests by point_in_multipolygon_rtree()
-*/
-typedef struct
-{
-	RTREE_NODE **ringIndices;
-	int* ringCounts;
-	int polyCount;
-} RTREE_POLY_CACHE;
-
-
-typedef struct
-{
-	GeomCache             gcache;
-	RTREE_POLY_CACHE      *index;
-} RTreeGeomCache;
-
-/**
-* Retrieves a collection of line segments given the root and crossing value.
-*/
-LWMLINE *RTreeFindLineSegments(RTREE_NODE *root, double value);
-
-
-/**
-* Checks for a cache hit against the provided geometry and returns
-* a pre-built index structure (RTREE_POLY_CACHE) if one exists. Otherwise
-* builds a new one and returns that.
-*/
-RTREE_POLY_CACHE *GetRtreeCache(FunctionCallInfo fcinfo, SHARED_GSERIALIZED *g1);
-
-int pip_short_circuit(RTREE_POLY_CACHE *poly_cache, LWPOINT *point, const GSERIALIZED *gpoly);
-
-#endif /* !defined _LWGEOM_RTREE_H */

-----------------------------------------------------------------------

Summary of changes:
 NEWS                                               |   1 +
 liblwgeom/Makefile.in                              |   2 +
 liblwgeom/cunit/cu_tree.c                          | 191 ++++++
 liblwgeom/intervaltree.c                           | 530 ++++++++++++++++
 liblwgeom/intervaltree.h                           |  76 +++
 liblwgeom/liblwgeom_internal.h                     |   1 +
 liblwgeom/lwutil.c                                 |  10 +-
 liblwgeom/ptarray.c                                |  12 +-
 libpgcommon/lwgeom_cache.c                         |   2 +-
 libpgcommon/lwgeom_cache.h                         |  24 +-
 postgis/Makefile.in                                |   9 +-
 postgis/lwgeom_functions_analytic.c                | 267 +-------
 postgis/lwgeom_geos.c                              | 242 ++------
 postgis/lwgeom_geos_prepared.h                     |  43 +-
 postgis/lwgeom_itree.c                             | 292 +++++++++
 ...{lwgeom_functions_analytic.h => lwgeom_itree.h} |  22 +-
 postgis/lwgeom_rectree.c                           |   2 +-
 postgis/lwgeom_rtree.c                             | 681 ---------------------
 postgis/lwgeom_rtree.h                             |  87 ---
 19 files changed, 1200 insertions(+), 1294 deletions(-)
 create mode 100644 liblwgeom/intervaltree.c
 create mode 100644 liblwgeom/intervaltree.h
 create mode 100644 postgis/lwgeom_itree.c
 rename postgis/{lwgeom_functions_analytic.h => lwgeom_itree.h} (59%)
 delete mode 100644 postgis/lwgeom_rtree.c
 delete mode 100644 postgis/lwgeom_rtree.h


hooks/post-receive
-- 
PostGIS