/**********************************************************************
 *
 * 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) 2009-2012 Paul Ramsey <pramsey@cleverelephant.ca>
 *
 **********************************************************************/

#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
#include "lwtree.h"
#include "measures.h"

static inline int
rect_node_is_leaf(const RECT_NODE *node)
{
	return node->type == RECT_NODE_LEAF_TYPE;
}

/*
* Support qsort of nodes for collection/multi types so nodes
* are in "spatial adjacent" order prior to merging.
*/
static int
rect_node_cmp(const void *pn1, const void *pn2)
{
	GBOX b1, b2;
	RECT_NODE *n1 = *((RECT_NODE**)pn1);
	RECT_NODE *n2 = *((RECT_NODE**)pn2);
	uint64_t h1, h2;
	b1.flags = 0;
	b1.xmin = n1->xmin;
	b1.xmax = n1->xmax;
	b1.ymin = n1->ymin;
	b1.ymax = n1->ymax;

	b2.flags = 0;
	b2.xmin = n2->xmin;
	b2.xmax = n2->xmax;
	b2.ymin = n2->ymin;
	b2.ymax = n2->ymax;

	h1 = gbox_get_sortable_hash(&b1, 0);
	h2 = gbox_get_sortable_hash(&b2, 0);
	return h1 < h2 ? -1 : (h1 > h2 ? 1 : 0);
}

/**
* Recurse from top of node tree and free all children.
* does not free underlying point array.
*/
void
rect_tree_free(RECT_NODE *node)
{
	int i;
	if (!node) return;
	if (!rect_node_is_leaf(node))
	{
		for (i = 0; i < node->i.num_nodes; i++)
		{
			rect_tree_free(node->i.nodes[i]);
			node->i.nodes[i] = NULL;
		}
	}
	lwfree(node);
}

static int
rect_leaf_node_intersects(RECT_NODE_LEAF *n1, RECT_NODE_LEAF *n2)
{
	const POINT2D *p1, *p2, *p3, *q1, *q2, *q3;
	DISTPTS dl;
	lw_dist2d_distpts_init(&dl, 1);
	switch (n1->seg_type)
	{
		case RECT_NODE_SEG_POINT:
		{
			p1 = getPoint2d_cp(n1->pa, n1->seg_num);

			switch (n2->seg_type)
			{
				case RECT_NODE_SEG_POINT:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					lw_dist2d_pt_pt(q1, p1, &dl);
					return dl.distance == 0.0;

				case RECT_NODE_SEG_LINEAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
					lw_dist2d_pt_seg(p1, q1, q2, &dl);
					return dl.distance == 0.0;

				case RECT_NODE_SEG_CIRCULAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
					q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
					lw_dist2d_pt_arc(p1, q1, q2, q3, &dl);
					return dl.distance == 0.0;

				default:
					lwerror("%s: unsupported segment type", __func__);
					break;
			}

			break;
		}

		case RECT_NODE_SEG_LINEAR:
		{
			p1 = getPoint2d_cp(n1->pa, n1->seg_num);
			p2 = getPoint2d_cp(n1->pa, n1->seg_num+1);

			switch (n2->seg_type)
			{
				case RECT_NODE_SEG_POINT:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					lw_dist2d_pt_seg(q1, p1, p2, &dl);
					return dl.distance == 0.0;

				case RECT_NODE_SEG_LINEAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
					return lw_segment_intersects(p1, p2, q1, q2) > 0;

				case RECT_NODE_SEG_CIRCULAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
					q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
					lw_dist2d_seg_arc(p1, p2, q1, q2, q3, &dl);
					return dl.distance == 0.0;

				default:
					lwerror("%s: unsupported segment type", __func__);
					break;
			}

			break;
		}
		case RECT_NODE_SEG_CIRCULAR:
		{
			p1 = getPoint2d_cp(n1->pa, n1->seg_num*2);
			p2 = getPoint2d_cp(n1->pa, n1->seg_num*2+1);
			p3 = getPoint2d_cp(n1->pa, n1->seg_num*2+2);

			switch (n2->seg_type)
			{
				case RECT_NODE_SEG_POINT:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					lw_dist2d_pt_arc(q1, p1, p2, p3, &dl);
					return dl.distance == 0.0;

				case RECT_NODE_SEG_LINEAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
					lw_dist2d_seg_arc(q1, q2, p1, p2, p3, &dl);
					return dl.distance == 0.0;

				case RECT_NODE_SEG_CIRCULAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
					q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
					lw_dist2d_arc_arc(p1, p2, p3, q1, q2, q3, &dl);
					return dl.distance == 0.0;

				default:
					lwerror("%s: unsupported segment type", __func__);
					break;
			}

			break;
		}
		default:
			return LW_FALSE;
	}
	return LW_FALSE;
}


/*
* Returns 1 if segment is to the right of point.
*/
static inline int
rect_leaf_node_segment_side(RECT_NODE_LEAF *node, const POINT2D *q, int *on_boundary)
{
	const POINT2D *p1, *p2, *p3;
	switch (node->seg_type)
	{
		case RECT_NODE_SEG_LINEAR:
		{
			int side;
			p1 = getPoint2d_cp(node->pa, node->seg_num);
			p2 = getPoint2d_cp(node->pa, node->seg_num+1);

			side = lw_segment_side(p1, p2, q);

			/* Always note case where we're on boundary */
			if (side == 0 && lw_pt_in_seg(q, p1, p2))
			{
				*on_boundary = LW_TRUE;
				return 0;
			}

			/* Segment points up and point is on left */
			if (p1->y < p2->y && side == -1 && q->y != p2->y)
			{
				return 1;
			}

			/* Segment points down and point is on right */
			if (p1->y > p2->y && side == 1 && q->y != p2->y)
			{
				return 1;
			}

			/* Segment is horizontal, do we cross first point? */
			if (p1->y == p2->y && q->x < p1->x)
			{
				return 1;
			}

			return 0;
		}
		case RECT_NODE_SEG_CIRCULAR:
		{
			int arc_side, seg_side;

			p1 = getPoint2d_cp(node->pa, node->seg_num*2);
			p2 = getPoint2d_cp(node->pa, node->seg_num*2+1);
			p3 = getPoint2d_cp(node->pa, node->seg_num*2+2);

			/* Always note case where we're on boundary */
			arc_side = lw_arc_side(p1, p2, p3, q);
			if (arc_side == 0)
			{
				*on_boundary = LW_TRUE;
				return 0;
			}

			seg_side = lw_segment_side(p1, p3, q);
			if (seg_side == arc_side)
			{
				/* Segment points up and point is on left */
				if (p1->y < p3->y && seg_side == -1 && q->y != p3->y)
				{
					return 1;
				}

				/* Segment points down and point is on right */
				if (p1->y > p3->y && seg_side == 1 && q->y != p3->y)
				{
					return 1;
				}
			}
			else
			{
				/* Segment points up and point is on left */
				if (p1->y < p3->y && seg_side == 1 && q->y != p3->y)
				{
					return 1;
				}

				/* Segment points down and point is on right */
				if (p1->y > p3->y && seg_side == -1 && q->y != p3->y)
				{
					return 1;
				}

				/* Segment is horizontal, do we cross first point? */
				if (p1->y == p3->y)
				{
					return 1;
				}
			}

			return 0;

		}
		default:
		{
			lwerror("%s: unsupported seg_type - %d", __func__, node->seg_type);
			return 0;
		}
	}

	return 0;
}

/*
* Only pass the head of a ring. Either a LinearRing from a polygon
* or a CompoundCurve from a CurvePolygon.
* Takes a horizontal line through the ring, and adds up the
* crossing directions. An "inside" point will be on the same side
* ("both left" or "both right") of the edges the line crosses,
* so it will have a non-zero return sum. An "outside" point will
* be on both sides, and will have a zero return sum.
*/
static int
rect_tree_ring_contains_point(RECT_NODE *node, const POINT2D *pt, int *on_boundary)
{
	/* Only test nodes that straddle our stabline vertically */
	/* and might be to the right horizontally */
	if (node->ymin <= pt->y && pt->y <= node->ymax && pt->x <= node->xmax)
	{
		if (rect_node_is_leaf(node))
		{
			return rect_leaf_node_segment_side(&node->l, pt, on_boundary);
		}
		else
		{
			int i, r = 0;
			for (i = 0; i < node->i.num_nodes; i++)
			{
				r += rect_tree_ring_contains_point(node->i.nodes[i], pt, on_boundary);
			}
			return r;
		}
	}
	return 0;
}

/*
* Only pass in the head of an "area" type. Polygon or CurvePolygon.
* Sums up containment of exterior (+1) and interior (-1) rings, so
* that zero is uncontained, +1 is contained and negative is an error
* (multiply contained by interior rings?)
*/
static int
rect_tree_area_contains_point(RECT_NODE *node, const POINT2D *pt)
{
	/* Can't do anything with a leaf node, makes no sense */
	if (rect_node_is_leaf(node))
		return 0;

	/* Iterate into area until we find ring heads */
	if (node->i.ring_type == RECT_NODE_RING_NONE)
	{
		int i, sum = 0;
		for (i = 0; i < node->i.num_nodes; i++)
			sum += rect_tree_area_contains_point(node->i.nodes[i], pt);
		return sum;
	}
	/* See if the ring encloses the point */
	else
	{
		int on_boundary = 0;
		int edge_crossing_count = rect_tree_ring_contains_point(node, pt, &on_boundary);
		/* Odd number of crossings => contained */
		int contained = (edge_crossing_count % 2 == 1);
		/* External rings return positive containment, interior ones negative, */
		/* so that a point-in-hole case nets out to zero (contained by both */
		/* interior and exterior rings. */
		if (node->i.ring_type == RECT_NODE_RING_INTERIOR)
		{
			return on_boundary ? 0 : -1 * contained;
		}
		else
		{
			return contained || on_boundary;
		}

	}
}

/*
* Simple containment test for node/point inputs
*/
static int
rect_node_bounds_point(RECT_NODE *node, const POINT2D *pt)
{
	if (pt->y < node->ymin || pt->y > node->ymax ||
		pt->x < node->xmin || pt->x > node->xmax)
		return LW_FALSE;
	else
		return LW_TRUE;
}

/*
* Pass in arbitrary tree, get back true if point is contained or on boundary,
* and false otherwise.
*/
int
rect_tree_contains_point(RECT_NODE *node, const POINT2D *pt)
{
	int i, c;

	/* Object cannot contain point if bounds don't */
	if (!rect_node_bounds_point(node, pt))
		return 0;

	switch (node->geom_type)
	{
		case POLYGONTYPE:
		case CURVEPOLYTYPE:
			return rect_tree_area_contains_point(node, pt) > 0;

		case MULTIPOLYGONTYPE:
		case MULTISURFACETYPE:
		case COLLECTIONTYPE:
		{
			for (i = 0; i < node->i.num_nodes; i++)
			{
				c = rect_tree_contains_point(node->i.nodes[i], pt);
				if (c) return LW_TRUE;
			}
			return LW_FALSE;
		}

		default:
			return LW_FALSE;
	}
}

/*
* For area types, doing intersects and distance, we will
* need to do a point-in-poly test first to find the full-contained
* case where an intersection exists without any edges actually
* intersecting.
*/
static int
rect_tree_is_area(const RECT_NODE *node)
{
	switch (node->geom_type)
	{
		case POLYGONTYPE:
		case CURVEPOLYTYPE:
		case MULTISURFACETYPE:
			return LW_TRUE;

		case COLLECTIONTYPE:
		{
			if (rect_node_is_leaf(node))
				return LW_FALSE;
			else
			{
				int i;
				for (i = 0; i < node->i.num_nodes; i++)
				{
					if (rect_tree_is_area(node->i.nodes[i]))
						return LW_TRUE;
				}
				return LW_FALSE;
			}
		}

		default:
			return LW_FALSE;
	}
}

static RECT_NODE_SEG_TYPE lwgeomTypeArc[] =
{
	RECT_NODE_SEG_UNKNOWN,  // "Unknown"
	RECT_NODE_SEG_POINT,    // "Point"
	RECT_NODE_SEG_LINEAR,   // "LineString"
	RECT_NODE_SEG_LINEAR,   // "Polygon"
	RECT_NODE_SEG_UNKNOWN,  // "MultiPoint"
	RECT_NODE_SEG_LINEAR,   // "MultiLineString"
	RECT_NODE_SEG_LINEAR,   // "MultiPolygon"
	RECT_NODE_SEG_UNKNOWN,  // "GeometryCollection"
	RECT_NODE_SEG_CIRCULAR, // "CircularString"
	RECT_NODE_SEG_UNKNOWN,  // "CompoundCurve"
	RECT_NODE_SEG_UNKNOWN,  // "CurvePolygon"
	RECT_NODE_SEG_UNKNOWN,  // "MultiCurve"
	RECT_NODE_SEG_UNKNOWN,  // "MultiSurface"
	RECT_NODE_SEG_LINEAR,   // "PolyhedralSurface"
	RECT_NODE_SEG_LINEAR,   // "Triangle"
	RECT_NODE_SEG_LINEAR    // "Tin"
};

/*
* Create a new leaf node.
*/
static RECT_NODE *
rect_node_leaf_new(const POINTARRAY *pa, int seg_num, int geom_type)
{
	const POINT2D *p1, *p2, *p3;
	RECT_NODE *node;
	GBOX gbox;
	RECT_NODE_SEG_TYPE seg_type = lwgeomTypeArc[geom_type];

	switch (seg_type)
	{
		case RECT_NODE_SEG_POINT:
		{
			p1 = getPoint2d_cp(pa, seg_num);
			gbox.xmin = gbox.xmax = p1->x;
			gbox.ymin = gbox.ymax = p1->y;
			break;
		}

		case RECT_NODE_SEG_LINEAR:
		{
			p1 = getPoint2d_cp(pa, seg_num);
			p2 = getPoint2d_cp(pa, seg_num+1);
			/* Zero length edge, doesn't get a node */
			if ((p1->x == p2->x) && (p1->y == p2->y))
				return NULL;
			gbox.xmin = FP_MIN(p1->x, p2->x);
			gbox.xmax = FP_MAX(p1->x, p2->x);
			gbox.ymin = FP_MIN(p1->y, p2->y);
			gbox.ymax = FP_MAX(p1->y, p2->y);
			break;
		}

		case RECT_NODE_SEG_CIRCULAR:
		{
			p1 = getPoint2d_cp(pa, 2*seg_num);
			p2 = getPoint2d_cp(pa, 2*seg_num+1);
			p3 = getPoint2d_cp(pa, 2*seg_num+2);
			/* Zero length edge, doesn't get a node */
			if ((p1->x == p2->x) && (p2->x == p3->x) &&
				(p1->y == p2->y) && (p2->y == p3->y))
				return NULL;
			lw_arc_calculate_gbox_cartesian_2d(p1, p2, p3, &gbox);
			break;
		}

		default:
		{
			lwerror("%s: unsupported seg_type - %d", __func__, seg_type);
			return NULL;
		}
	}

	node = lwalloc(sizeof(RECT_NODE));
	node->type = RECT_NODE_LEAF_TYPE;
	node->geom_type = geom_type;
	node->xmin = gbox.xmin;
	node->xmax = gbox.xmax;
	node->ymin = gbox.ymin;
	node->ymax = gbox.ymax;
	node->l.seg_num = seg_num;
	node->l.seg_type = seg_type;
	node->l.pa = pa;
	return node;
}


static void
rect_node_internal_add_node(RECT_NODE *node, RECT_NODE *add)
{
	if (rect_node_is_leaf(node))
		lwerror("%s: call on leaf node", __func__);
	node->xmin = FP_MIN(node->xmin, add->xmin);
	node->xmax = FP_MAX(node->xmax, add->xmax);
	node->ymin = FP_MIN(node->ymin, add->ymin);
	node->ymax = FP_MAX(node->ymax, add->ymax);
	node->i.nodes[node->i.num_nodes++] = add;
	return;
}


static RECT_NODE *
rect_node_internal_new(const RECT_NODE *seed)
{
	RECT_NODE *node = lwalloc(sizeof(RECT_NODE));
	node->xmin = seed->xmin;
	node->xmax = seed->xmax;
	node->ymin = seed->ymin;
	node->ymax = seed->ymax;
	node->geom_type = seed->geom_type;
	node->type = RECT_NODE_INTERNAL_TYPE;
	node->i.num_nodes = 0;
	node->i.ring_type = RECT_NODE_RING_NONE;
	node->i.sorted = 0;
	return node;
}

/*
* We expect the incoming nodes to be in a spatially coherent
* order. For incoming nodes derived from point arrays,
* the very fact that they are
* a vertex list implies a reasonable ordering: points nearby in
* the list will be nearby in space. For incoming nodes from higher
* level structures (collections, etc) the caller should sort the
* nodes using a z-order first, so that this merge step results in a
* spatially coherent structure.
*/
static RECT_NODE *
rect_nodes_merge(RECT_NODE ** nodes, uint32_t num_nodes)
{
	if (num_nodes < 1)
	{
		return NULL;
	}

	while (num_nodes > 1)
	{
		uint32_t i, k = 0;
		RECT_NODE *node = NULL;
		for (i = 0; i < num_nodes; i++)
		{
			if (!node)
				node = rect_node_internal_new(nodes[i]);

			rect_node_internal_add_node(node, nodes[i]);

			if (node->i.num_nodes == RECT_NODE_SIZE)
			{
				nodes[k++] = node;
				node = NULL;
			}
		}
		if (node)
			nodes[k++] = node;
		num_nodes = k;
	}

	return nodes[0];
}

/*
* Build a tree of nodes from a point array, one node per edge.
*/
RECT_NODE *
rect_tree_from_ptarray(const POINTARRAY *pa, int geom_type)
{
	int num_edges = 0, i = 0, j = 0;
	RECT_NODE_SEG_TYPE seg_type = lwgeomTypeArc[geom_type];
	RECT_NODE **nodes = NULL;
	RECT_NODE *tree = NULL;

	/* No-op on empty ring/line/pt */
	if ( pa->npoints < 1 )
		return NULL;

	/* For arcs, 3 points per edge, for lines, 2 per edge */
	switch(seg_type)
	{
		case RECT_NODE_SEG_POINT:
			return rect_node_leaf_new(pa, 0, geom_type);
			break;
		case RECT_NODE_SEG_LINEAR:
			num_edges = pa->npoints - 1;
			break;
		case RECT_NODE_SEG_CIRCULAR:
			num_edges = (pa->npoints - 1)/2;
			break;
		default:
			lwerror("%s: unsupported seg_type - %d", __func__, seg_type);
	}

	/* First create a flat list of nodes, one per edge. */
	nodes = lwalloc(sizeof(RECT_NODE*) * num_edges);
	for (i = 0; i < num_edges; i++)
	{
		RECT_NODE *node = rect_node_leaf_new(pa, i, geom_type);
		if (node) /* Not zero length? */
			nodes[j++] = node;
	}

	/* Merge the list into a tree */
	tree = rect_nodes_merge(nodes, j);

	/* Free the old list structure, leaving the tree in place */
	lwfree(nodes);

	/* Return top of tree */
	return tree;
}

LWGEOM *
rect_tree_to_lwgeom(const RECT_NODE *node)
{
	LWGEOM *poly = (LWGEOM*)lwpoly_construct_envelope(0, node->xmin, node->ymin, node->xmax, node->ymax);
	if (rect_node_is_leaf(node))
	{
		return poly;
	}
	else
	{
		int i;
		LWCOLLECTION *col = lwcollection_construct_empty(COLLECTIONTYPE, 0, 0, 0);
		lwcollection_add_lwgeom(col, poly);
		for (i = 0; i < node->i.num_nodes; i++)
		{
			lwcollection_add_lwgeom(col, rect_tree_to_lwgeom(node->i.nodes[i]));
		}
		return (LWGEOM*)col;
	}
}

char *
rect_tree_to_wkt(const RECT_NODE *node)
{
	LWGEOM *geom = rect_tree_to_lwgeom(node);
	char *wkt = lwgeom_to_wkt(geom, WKT_ISO, 12, 0);
	lwgeom_free(geom);
	return wkt;
}

void
rect_tree_printf(const RECT_NODE *node, int depth)
{
	printf("%*s----\n", depth, "");
	printf("%*stype: %d\n", depth, "", node->type);
	printf("%*sgeom_type: %d\n", depth, "", node->geom_type);
	printf("%*sbox: %g %g, %g %g\n", depth, "", node->xmin, node->ymin, node->xmax, node->ymax);
	if (node->type == RECT_NODE_LEAF_TYPE)
	{
		printf("%*sseg_type: %d\n", depth, "", node->l.seg_type);
		printf("%*sseg_num: %d\n", depth, "", node->l.seg_num);
	}
	else
	{
		int i;
		for (i = 0; i < node->i.num_nodes; i++)
		{
			rect_tree_printf(node->i.nodes[i], depth+2);
		}
	}
}

static RECT_NODE *
rect_tree_from_lwpoint(const LWGEOM *lwgeom)
{
	const LWPOINT *lwpt = (const LWPOINT*)lwgeom;
	return rect_tree_from_ptarray(lwpt->point, lwgeom->type);
}

static RECT_NODE *
rect_tree_from_lwline(const LWGEOM *lwgeom)
{
	const LWLINE *lwline = (const LWLINE*)lwgeom;
	return rect_tree_from_ptarray(lwline->points, lwgeom->type);
}

static RECT_NODE *
rect_tree_from_lwpoly(const LWGEOM *lwgeom)
{
	RECT_NODE **nodes;
	RECT_NODE *tree;
	uint32_t i, j = 0;
	const LWPOLY *lwpoly = (const LWPOLY*)lwgeom;

	if (lwpoly->nrings < 1)
		return NULL;

	nodes = lwalloc(sizeof(RECT_NODE*) * lwpoly->nrings);
	for (i = 0; i < lwpoly->nrings; i++)
	{
		RECT_NODE *node = rect_tree_from_ptarray(lwpoly->rings[i], lwgeom->type);
		if (node)
		{
			node->i.ring_type = i ? RECT_NODE_RING_INTERIOR : RECT_NODE_RING_EXTERIOR;
			nodes[j++] = node;
		}
	}
	tree = rect_nodes_merge(nodes, j);
	tree->geom_type = lwgeom->type;
	lwfree(nodes);
	return tree;
}

static RECT_NODE *
rect_tree_from_lwcurvepoly(const LWGEOM *lwgeom)
{
	RECT_NODE **nodes;
	RECT_NODE *tree;
	uint32_t i, j = 0;
	const LWCURVEPOLY *lwcol = (const LWCURVEPOLY*)lwgeom;

	if (lwcol->nrings < 1)
		return NULL;

	nodes = lwalloc(sizeof(RECT_NODE*) * lwcol->nrings);
	for (i = 0; i < lwcol->nrings; i++)
	{
		RECT_NODE *node = rect_tree_from_lwgeom(lwcol->rings[i]);
		if (node)
		{
			/*
			* In the case of arc circle, it's possible for a ring to consist
			* of a single closed edge. That will arrive as a leaf node. We
			* need to wrap that node in an internal node with an appropriate
			* ring type so all the other code can try and make sense of it.
			*/
			if (node->type == RECT_NODE_LEAF_TYPE)
			{
				RECT_NODE *internal = rect_node_internal_new(node);
				rect_node_internal_add_node(internal, node);
				node = internal;
			}
			/* Each subcomponent is a ring */
			node->i.ring_type = i ? RECT_NODE_RING_INTERIOR : RECT_NODE_RING_EXTERIOR;
			nodes[j++] = node;
		}
	}
	/* Put the top nodes in a z-order curve for a spatially coherent */
	/* tree after node merge */
	qsort(nodes, j, sizeof(RECT_NODE*), rect_node_cmp);

	tree = rect_nodes_merge(nodes, j);

	tree->geom_type = lwgeom->type;
	lwfree(nodes);
	return tree;

}

static RECT_NODE *
rect_tree_from_lwcollection(const LWGEOM *lwgeom)
{
	RECT_NODE **nodes;
	RECT_NODE *tree;
	uint32_t i, j = 0;
	const LWCOLLECTION *lwcol = (const LWCOLLECTION*)lwgeom;

	if (lwcol->ngeoms < 1)
		return NULL;

	/* Build one tree for each sub-geometry, then below */
	/* we merge the root notes of those trees to get a single */
	/* top node for the collection */
	nodes = lwalloc(sizeof(RECT_NODE*) * lwcol->ngeoms);
	for (i = 0; i < lwcol->ngeoms; i++)
	{
		RECT_NODE *node = rect_tree_from_lwgeom(lwcol->geoms[i]);
		if (node)
		{
			/* Curvepolygons are collections where the sub-geometries */
			/* are the rings, and will need to doint point-in-poly */
			/* tests in order to do intersects and distance calculations */
			/* correctly */
			if (lwgeom->type == CURVEPOLYTYPE)
				node->i.ring_type = i ? RECT_NODE_RING_INTERIOR : RECT_NODE_RING_EXTERIOR;
			nodes[j++] = node;
		}
	}
	/* Sort the nodes using a z-order curve, so that merging the nodes */
	/* gives a spatially coherent tree (near things are in near nodes) */
	/* Note: CompoundCurve has edges already spatially organized, no */
	/* sorting needed */
	if (lwgeom->type != COMPOUNDTYPE)
		qsort(nodes, j, sizeof(RECT_NODE*), rect_node_cmp);

	tree = rect_nodes_merge(nodes, j);

	tree->geom_type = lwgeom->type;
	lwfree(nodes);
	return tree;
}

RECT_NODE *
rect_tree_from_lwgeom(const LWGEOM *lwgeom)
{
	switch(lwgeom->type)
	{
		case POINTTYPE:
			return rect_tree_from_lwpoint(lwgeom);
		case TRIANGLETYPE:
		case CIRCSTRINGTYPE:
		case LINETYPE:
			return rect_tree_from_lwline(lwgeom);
		case POLYGONTYPE:
			return rect_tree_from_lwpoly(lwgeom);
		case CURVEPOLYTYPE:
			return rect_tree_from_lwcurvepoly(lwgeom);
		case COMPOUNDTYPE:
		case MULTICURVETYPE:
		case MULTISURFACETYPE:
		case MULTIPOINTTYPE:
		case MULTILINETYPE:
		case MULTIPOLYGONTYPE:
		case POLYHEDRALSURFACETYPE:
		case TINTYPE:
		case COLLECTIONTYPE:
			return rect_tree_from_lwcollection(lwgeom);
		default:
			lwerror("%s: Unknown geometry type: %s", __func__, lwtype_name(lwgeom->type));
			return NULL;
	}
	return NULL;
}

/*
* Get an actual coordinate point from a tree to use
* for point-in-polygon testing.
*/
static const POINT2D *
rect_tree_get_point(const RECT_NODE *node)
{
	if (!node) return NULL;
	if (rect_node_is_leaf(node))
		return getPoint2d_cp(node->l.pa, 0);
	else
		return rect_tree_get_point(node->i.nodes[0]);
}

static inline int
rect_node_intersects(const RECT_NODE *n1, const RECT_NODE *n2)
{
	if (n1->xmin > n2->xmax || n2->xmin > n1->xmax ||
		n1->ymin > n2->ymax || n2->ymin > n1->ymax)
	{
		return 0;
	}
	else
	{
		return 1;
	}
}

#if POSTGIS_DEBUG_LEVEL >= 4
static char *
rect_node_to_str(const RECT_NODE *n)
{
	char *buf = lwalloc(256);
	snprintf(buf, 256, "(%.9g %.9g,%.9g %.9g)",
		n->xmin, n->ymin, n->xmax, n->ymax);
   	return buf;
}
#endif

/*
* Work down to leaf nodes, until we find a pair of leaf nodes
* that intersect. Prune branches that do not intersect.
*/
static int
rect_tree_intersects_tree_recursive(RECT_NODE *n1, RECT_NODE *n2)
{
	int i, j;
#if POSTGIS_DEBUG_LEVEL >= 4
	char *n1_str = rect_node_to_str(n1);
	char *n2_str = rect_node_to_str(n2);
	LWDEBUGF(4,"n1 %s  n2 %s", n1, n2);
	lwfree(n1_str);
	lwfree(n2_str);
#endif
	/* There can only be an edge intersection if the rectangles overlap */
	if (rect_node_intersects(n1, n2))
	{
		LWDEBUG(4," interaction found");
		/* We can only test for a true intersection if the nodes are both leaf nodes */
		if (rect_node_is_leaf(n1) && rect_node_is_leaf(n2))
		{
			LWDEBUG(4,"  leaf node test");
			/* Check for true intersection */
			return rect_leaf_node_intersects(&n1->l, &n2->l);
		}
		else if (rect_node_is_leaf(n2) && !rect_node_is_leaf(n1))
		{
			for (i = 0; i < n1->i.num_nodes; i++)
			{
				if (rect_tree_intersects_tree_recursive(n1->i.nodes[i], n2))
					return LW_TRUE;
			}
		}
		else if (rect_node_is_leaf(n1) && !rect_node_is_leaf(n2))
		{
			for (i = 0; i < n2->i.num_nodes; i++)
			{
				if (rect_tree_intersects_tree_recursive(n2->i.nodes[i], n1))
					return LW_TRUE;
			}
		}
		else
		{
			for (j = 0; j < n1->i.num_nodes; j++)
			{
				for (i = 0; i < n2->i.num_nodes; i++)
				{
					if (rect_tree_intersects_tree_recursive(n2->i.nodes[i], n1->i.nodes[j]))
						return LW_TRUE;
				}
			}
		}
	}
	return LW_FALSE;
}


int
rect_tree_intersects_tree(RECT_NODE *n1, RECT_NODE *n2)
{
	/*
	* It is possible for an area to intersect another object
	* without any edges intersecting, if the object is fully contained.
	* If that is so, then any point in the object will be contained,
	* so we do a quick point-in-poly test first for those cases
	*/
	if (rect_tree_is_area(n1) &&
		rect_tree_contains_point(n1, rect_tree_get_point(n2)))
	{
		return LW_TRUE;
	}

	if (rect_tree_is_area(n2) &&
		rect_tree_contains_point(n2, rect_tree_get_point(n1)))
	{
		return LW_TRUE;
	}

	/*
	* Not contained, so intersection can only happen if
	* edges actually intersect.
	*/
	return rect_tree_intersects_tree_recursive(n1, n2);
}

/*
* For slightly more speed when doing distance comparisons,
* where we only care if d1 < d2 and not what the actual value
* of d1 or d2 is, we use the squared distance and avoid the
* sqrt()
*/
static inline double
distance_sq(double x1, double y1, double x2, double y2)
{
	double dx = x1-x2;
	double dy = y1-y2;
	return dx*dx + dy*dy;
}

static inline double
distance(double x1, double y1, double x2, double y2)
{
	return sqrt(distance_sq(x1, y1, x2, y2));
}

/*
* The closest any two objects in two nodes can be is the smallest
* distance between the nodes themselves.
*/
static inline double
rect_node_min_distance(const RECT_NODE *n1, const RECT_NODE *n2)
{
	int   left = n1->xmin > n2->xmax;
	int  right = n1->xmax < n2->xmin;
	int bottom = n1->ymin > n2->ymax;
	int    top = n1->ymax < n2->ymin;

	//lwnotice("rect_node_min_distance");

	if (top && left)
		return distance(n1->xmin, n1->ymax, n2->xmax, n2->ymin);
	else if (top && right)
		return distance(n1->xmax, n1->ymax, n2->xmin, n2->ymin);
	else if (bottom && left)
		return distance(n1->xmin, n1->ymin, n2->xmax, n2->ymax);
	else if (bottom && right)
		return distance(n1->xmax, n1->ymin, n2->xmin, n2->ymax);
	else if (left)
		return n1->xmin - n2->xmax;
	else if (right)
		return n2->xmin - n1->xmax;
	else if (bottom)
		return n1->ymin - n2->ymax;
	else if (top)
		return n2->ymin - n1->ymax;
	else
		return 0.0;

	return 0.0;
}

/*
* The furthest apart the objects in two nodes can be is if they
* are at opposite corners of the bbox that contains both nodes
*/
static inline double
rect_node_max_distance(const RECT_NODE *n1, const RECT_NODE *n2)
{
	double xmin = FP_MIN(n1->xmin, n2->xmin);
	double ymin = FP_MIN(n1->ymin, n2->ymin);
	double xmax = FP_MAX(n1->xmax, n2->xmax);
	double ymax = FP_MAX(n1->ymax, n2->ymax);
	double dx = xmax - xmin;
	double dy = ymax - ymin;
	//lwnotice("rect_node_max_distance");
	return sqrt(dx*dx + dy*dy);
}

/*
* Leaf nodes represent individual edges from the original shape.
* As such, they can be either points (if original was a (multi)point)
* two-point straight edges (for linear types), or
* three-point circular arcs (for curvilinear types).
* The distance calculations for each possible combination of primitive
* edges is different, so there's a big case switch in here to match
* up the right combination of inputs to the right distance calculation.
*/
static double
rect_leaf_node_distance(const RECT_NODE_LEAF *n1, const RECT_NODE_LEAF *n2, RECT_TREE_DISTANCE_STATE *state)
{
	const POINT2D *p1, *p2, *p3, *q1, *q2, *q3;
	DISTPTS dl;

	//lwnotice("rect_leaf_node_distance, %d<->%d", n1->seg_num, n2->seg_num);

	lw_dist2d_distpts_init(&dl, DIST_MIN);

	switch (n1->seg_type)
	{
		case RECT_NODE_SEG_POINT:
		{
			p1 = getPoint2d_cp(n1->pa, n1->seg_num);

			switch (n2->seg_type)
			{
				case RECT_NODE_SEG_POINT:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					lw_dist2d_pt_pt(q1, p1, &dl);
					break;

				case RECT_NODE_SEG_LINEAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
					lw_dist2d_pt_seg(p1, q1, q2, &dl);
					break;

				case RECT_NODE_SEG_CIRCULAR:
				q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
					q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
					lw_dist2d_pt_arc(p1, q1, q2, q3, &dl);
					break;

				default:
					lwerror("%s: unsupported segment type", __func__);
			}
			break;
		}

		case RECT_NODE_SEG_LINEAR:
		{
			p1 = getPoint2d_cp(n1->pa, n1->seg_num);
			p2 = getPoint2d_cp(n1->pa, n1->seg_num+1);

			switch (n2->seg_type)
			{
				case RECT_NODE_SEG_POINT:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					lw_dist2d_pt_seg(q1, p1, p2, &dl);
					break;

				case RECT_NODE_SEG_LINEAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
					lw_dist2d_seg_seg(q1, q2, p1, p2, &dl);
					// lwnotice(
					// 	"%d\tLINESTRING(%g %g,%g %g)\t%d\tLINESTRING(%g %g,%g %g)\t%g\t%g\t%g",
					// 	n1->seg_num,
					// 	p1->x, p1->y, p2->x, p2->y,
					// 	n2->seg_num,
					// 	q1->x, q1->y, q2->x, q2->y,
					// 	dl.distance, state->min_dist, state->max_dist);
					break;

				case RECT_NODE_SEG_CIRCULAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
					q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
					lw_dist2d_seg_arc(p1, p2, q1, q2, q3, &dl);
					break;

				default:
					lwerror("%s: unsupported segment type", __func__);
			}
			break;
		}
		case RECT_NODE_SEG_CIRCULAR:
		{
			p1 = getPoint2d_cp(n1->pa, n1->seg_num*2);
			p2 = getPoint2d_cp(n1->pa, n1->seg_num*2+1);
			p3 = getPoint2d_cp(n1->pa, n1->seg_num*2+2);

			switch (n2->seg_type)
			{
				case RECT_NODE_SEG_POINT:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					lw_dist2d_pt_arc(q1, p1, p2, p3, &dl);
					break;

				case RECT_NODE_SEG_LINEAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
					lw_dist2d_seg_arc(q1, q2, p1, p2, p3, &dl);
					break;

				case RECT_NODE_SEG_CIRCULAR:
					q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
					q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
					q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
					lw_dist2d_arc_arc(p1, p2, p3, q1, q2, q3, &dl);
					break;

				default:
					lwerror("%s: unsupported segment type", __func__);
			}
			break;
		}
		default:
			lwerror("%s: unsupported segment type", __func__);
	}

	/* If this is a new global minima, save it */
	if (dl.distance < state->min_dist)
	{
		state->min_dist = dl.distance;
		state->p1 = dl.p1;
		state->p2 = dl.p2;
	}

	return dl.distance;
}


static int
rect_tree_node_sort_cmp(const void *a, const void *b)
{
	RECT_NODE *n1 = *(RECT_NODE**)a;
	RECT_NODE *n2 = *(RECT_NODE**)b;
	if (n1->d < n2->d) return -1;
	else if (n1->d > n2->d) return 1;
	else return 0;
}

/* Calculate the center of a node */
static inline void
rect_node_to_p2d(const RECT_NODE *n, POINT2D *pt)
{
	pt->x = (n->xmin + n->xmax)/2;
	pt->y = (n->ymin + n->ymax)/2;
}

/*
* (If necessary), sort the children of each node in
* order of their distance from the enter of the other node.
*/
static void
rect_tree_node_sort(RECT_NODE *n1, RECT_NODE *n2)
{
	int i;
	RECT_NODE *n;
	POINT2D c1, c2, c;
	if (!rect_node_is_leaf(n1) && ! n1->i.sorted)
	{
		rect_node_to_p2d(n2, &c2);
		/* Distance of each child from center of other node */
		for (i = 0; i < n1->i.num_nodes; i++)
		{
			n = n1->i.nodes[i];
			rect_node_to_p2d(n, &c);
			n->d = distance2d_sqr_pt_pt(&c2, &c);
		}
		/* Sort the children by distance */
		n1->i.sorted = 1;
		qsort(n1->i.nodes,
		         n1->i.num_nodes,
		         sizeof(RECT_NODE*),
		         rect_tree_node_sort_cmp);
	}
	if (!rect_node_is_leaf(n2) && ! n2->i.sorted)
	{
		rect_node_to_p2d(n1, &c1);
		/* Distance of each child from center of other node */
		for (i = 0; i < n2->i.num_nodes; i++)
		{
			n = n2->i.nodes[i];
			rect_node_to_p2d(n, &c);
			n->d = distance2d_sqr_pt_pt(&c1, &c);
		}
		/* Sort the children by distance */
		n2->i.sorted = 1;
		qsort(n2->i.nodes,
		         n2->i.num_nodes,
		         sizeof(RECT_NODE*),
		         rect_tree_node_sort_cmp);
	}
}

static double
rect_tree_distance_tree_recursive(RECT_NODE *n1, RECT_NODE *n2, RECT_TREE_DISTANCE_STATE *state)
{
	double min, max;

	/* Short circuit if we've already hit the minimum */
	if (state->min_dist < state->threshold || state->min_dist == 0.0)
		return state->min_dist;

	/* If your minimum is greater than anyone's maximum, you can't hold the winner */
	min = rect_node_min_distance(n1, n2);
	if (min > state->max_dist)
	{
		//lwnotice("pruning pair %p, %p", n1, n2);
		LWDEBUGF(4, "pruning pair %p, %p", n1, n2);
		return FLT_MAX;
	}

	/* If your maximum is a new low, we'll use that as our new global tolerance */
	max = rect_node_max_distance(n1, n2);
	if (max < state->max_dist)
		state->max_dist = max;

	/* Both leaf nodes, do a real distance calculation */
	if (rect_node_is_leaf(n1) && rect_node_is_leaf(n2))
	{
		return rect_leaf_node_distance(&n1->l, &n2->l, state);
	}
	/* Recurse into nodes */
	else
	{
		int i, j;
		double d_min = FLT_MAX;
		rect_tree_node_sort(n1, n2);
		if (rect_node_is_leaf(n1) && !rect_node_is_leaf(n2))
		{
			for (i = 0; i < n2->i.num_nodes; i++)
			{
				min = rect_tree_distance_tree_recursive(n1, n2->i.nodes[i], state);
				d_min = FP_MIN(d_min, min);
			}
		}
		else if (rect_node_is_leaf(n2) && !rect_node_is_leaf(n1))
		{
			for (i = 0; i < n1->i.num_nodes; i++)
			{
				min = rect_tree_distance_tree_recursive(n1->i.nodes[i], n2, state);
				d_min = FP_MIN(d_min, min);
			}
		}
		else
		{
			for (i = 0; i < n1->i.num_nodes; i++)
			{
				for (j = 0; j < n2->i.num_nodes; j++)
				{
					min = rect_tree_distance_tree_recursive(n1->i.nodes[i], n2->i.nodes[j], state);
					d_min = FP_MIN(d_min, min);
				}
			}
		}
		return d_min;
	}
}

double rect_tree_distance_tree(RECT_NODE *n1, RECT_NODE *n2, double threshold)
{
	double distance;
	RECT_TREE_DISTANCE_STATE state;

	/*
	* It is possible for an area to intersect another object
	* without any edges intersecting, if the object is fully contained.
	* If that is so, then any point in the object will be contained,
	* so we do a quick point-in-poly test first for those cases
	*/
	if (rect_tree_is_area(n1) &&
		rect_tree_contains_point(n1, rect_tree_get_point(n2)))
	{
		return 0.0;
	}

	if (rect_tree_is_area(n2) &&
		rect_tree_contains_point(n2, rect_tree_get_point(n1)))
	{
		return 0.0;
	}

	state.threshold = threshold;
	state.min_dist = FLT_MAX;
	state.max_dist = FLT_MAX;
	distance = rect_tree_distance_tree_recursive(n1, n2, &state);
	// *p1 = state.p1;
	// *p2 = state.p2;
	return distance;
}