Add an internal geometry tree for use in a native prepared geometry scheme for... (a6c08ab1) · Commits · postgis / postgis

liblwgeom/Makefile.in

+2 −1

Original line number	Diff line number	Diff line
		@@ -50,7 +50,8 @@ SA_OBJS = \
		g_serialized.o \
		g_util.o \
		lwgeodetic.o \
		lwspheroid.o
		lwspheroid.o \
		lwtree.o

		SA_HEADERS = \
		liblwgeom.h \

liblwgeom/lwalgorithm.c

+22 −19

Original line number	Diff line number	Diff line
		@@ -21,26 +21,30 @@
		** Return > 0.0 if point Q is right of segment P
		** Return = 0.0 if point Q in on segment P
		*/
		double lw_segment_side(POINT2D p1, POINT2D p2, POINT2D q)
		double lw_segment_side(const POINT2D p1, const POINT2D p2, const POINT2D *q)
		{
		return ( (q.x - p1.x) * (p2.y - p1.y) - (p2.x - p1.x) * (q.y - p1.y) );
		double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
		if( FP_IS_ZERO(side) )
		return 0.0;
		else
		return side;
		}


		int lw_segment_envelope_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
		int lw_segment_envelope_intersects(const POINT2D p1, const POINT2D p2, const POINT2D q1, const POINT2D q2)
		{
		double minq=LW_MIN(q1.x,q2.x);
		double maxq=LW_MAX(q1.x,q2.x);
		double minp=LW_MIN(p1.x,p2.x);
		double maxp=LW_MAX(p1.x,p2.x);
		double minq=LW_MIN(q1->x,q2->x);
		double maxq=LW_MAX(q1->x,q2->x);
		double minp=LW_MIN(p1->x,p2->x);
		double maxp=LW_MAX(p1->x,p2->x);

		if (FP_GT(minp,maxq) \|\| FP_LT(maxp,minq))
		return LW_FALSE;

		minq=LW_MIN(q1.y,q2.y);
		maxq=LW_MAX(q1.y,q2.y);
		minp=LW_MIN(p1.y,p2.y);
		maxp=LW_MAX(p1.y,p2.y);
		minq=LW_MIN(q1->y,q2->y);
		maxq=LW_MAX(q1->y,q2->y);
		minp=LW_MIN(p1->y,p2->y);
		maxp=LW_MAX(p1->y,p2->y);

		if (FP_GT(minp,maxq) \|\| FP_LT(maxp,minq))
		return LW_FALSE;
		@@ -62,7 +66,7 @@ int lw_segment_envelope_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q
		** SEG_CROSS_LEFT = 2,
		** SEG_CROSS_RIGHT = 3,
		*/
		int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
		int lw_segment_intersects(const POINT2D p1, const POINT2D p2, const POINT2D q1, const POINT2D q2)
		{

		double pq1, pq2, qp1, qp2;
		@@ -84,13 +88,13 @@ int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
		/* Are the start and end points of p on the same side of q? */
		qp1=lw_segment_side(q1,q2,p1);
		qp2=lw_segment_side(q1,q2,p2);
		if ((qp1>0 && qp2>0) \|\| (qp1<0 && qp2<0))
		if ( (qp1 > 0.0 && qp2 > 0.0) \|\| (qp1 < 0.0 && qp2 < 0.0) )
		{
		return SEG_NO_INTERSECTION;
		}

		/* Nobody is on one side or another? Must be colinear. */
		if (FP_IS_ZERO(pq1) && FP_IS_ZERO(pq2) && FP_IS_ZERO(qp1) && FP_IS_ZERO(qp2))
		if ( pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0 )
		{
		return SEG_COLINEAR;
		}
		@@ -104,13 +108,13 @@ int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
		LWDEBUGF(4, "qp1=%.15g qp2=%.15g", qp1, qp2);

		/* Second point of p or q touches, it's not a crossing. */
		if ( FP_IS_ZERO(pq2) \|\| FP_IS_ZERO(qp2) )
		if ( pq2 == 0.0 \|\| qp2 == 0.0 )
		{
		return SEG_NO_INTERSECTION;
		}

		/* First point of p touches, it's a "crossing". */
		if ( FP_IS_ZERO(pq1) )
		if ( pq1 == 0.0 )
		{
		if ( FP_GT(pq2,0.0) )
		return SEG_CROSS_RIGHT;
		@@ -119,7 +123,7 @@ int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
		}

		/* First point of q touches, it's a crossing. */
		if ( FP_IS_ZERO(qp1) )
		if ( qp1 == 0.0 )
		{
		if ( FP_LT(pq1,pq2) )
		return SEG_CROSS_RIGHT;
		@@ -135,7 +139,6 @@ int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)

		/* This should never happen! */
		return SEG_ERROR;

		}

		/**
		@@ -190,7 +193,7 @@ int lwline_crossing_direction(LWLINE l1, LWLINE l2)
		/* Update second point of p to next value */
		rv = getPoint2d_p(pa1, j, &p2);

		this_cross = lw_segment_intersects(p1, p2, q1, q2);
		this_cross = lw_segment_intersects(&p1, &p2, &q1, &q2);

		LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", this_cross, i, j, p1.x, p1.y, p2.x, p2.y);

liblwgeom/lwalgorithm.h

+3 −3

Original line number	Diff line number	Diff line
		@@ -23,9 +23,9 @@ enum CG_SEGMENT_INTERSECTION_TYPE {
		SEG_TOUCH_RIGHT = 5
		};

		double lw_segment_side(POINT2D p1, POINT2D p2, POINT2D q);
		int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2);
		int lw_segment_envelope_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2);
		double lw_segment_side(const POINT2D p1, const POINT2D p2, const POINT2D *q);
		int lw_segment_intersects(const POINT2D p1, const POINT2D p2, const POINT2D q1, const POINT2D q2);
		int lw_segment_envelope_intersects(const POINT2D p1, const POINT2D p2, const POINT2D q1, const POINT2D q2);


		enum CG_LINE_CROSS_TYPE {

liblwgeom/lwgeodetic.h

+9 −0

Original line number	Diff line number	Diff line
		@@ -34,6 +34,15 @@ typedef struct
		GEOGRAPHIC_POINT end;
		} GEOGRAPHIC_EDGE;

		/**
		* Holder for sorting points in distance algorithm
		*/
		typedef struct
		{
		double measure;
		uint32 index;
		} DISTANCE_ORDER;

		/**
		* Conversion functions
		*/

liblwgeom/lwtree.c

0 → 100644

+218 −0

Original line number	Diff line number	Diff line
		#include "lwalgorithm.h"
		#include "lwtree.h"


		/**
		* Internal nodes have their point references set to NULL.
		*/
		static int rect_node_is_leaf(const RECT_NODE *node)
		{
		return (node->p1 != NULL);
		}

		/**
		* Recurse from top of node tree and free all children.
		* does not free underlying point array.
		*/
		void rect_tree_free(RECT_NODE *node)
		{
		if( node->left_node )
		{
		rect_tree_free(node->left_node);
		node->left_node = 0;
		}
		if( node->right_node )
		{
		rect_tree_free(node->right_node);
		node->right_node = 0;
		}
		lwfree(node);
		}

		/* 0 => no containment */
		int rect_tree_contains_point(const RECT_NODE node, const POINT2D pt, int *on_boundary)
		{
		if( FP_CONTAINS_INCL(node->ymin, pt->y, node->ymax) )
		{
		if( rect_node_is_leaf(node) )
		{
		double side = lw_segment_side(node->p1, node->p2, pt);
		if ( side == 0.0 )
		*on_boundary = LW_TRUE;
		return (side < 0.0 ? -1 : 1 );
		}
		else
		{
		return rect_tree_contains_point(node->left_node, pt, on_boundary) +
		rect_tree_contains_point(node->right_node, pt, on_boundary);
		}
		}
		//printf("NOT in measure range\n");
		return 0;
		}

		int rect_tree_intersects_tree(const RECT_NODE n1, const RECT_NODE n2)
		{
		/* There can only be an edge intersection if the rectangles overlap */
		if( ! ( FP_GTEQ(n1->xmin, n2->xmax) \|\| FP_GTEQ(n2->xmin, n1->xmax) \|\| FP_GTEQ(n1->ymin, n2->ymax) \|\| FP_GTEQ(n2->ymin, n1->ymax) ) )
		{
		/* 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) )
		{
		/* Check for true intersection */
		if( lw_segment_intersects(n1->p1, n1->p2, n2->p1, n2->p2) )
		return LW_TRUE;
		else
		return LW_FALSE;
		}
		else
		{
		/* Recurse to children */
		if( rect_node_is_leaf(n1) )
		{
		if( rect_tree_intersects_tree(n2->left_node, n1) \|\| rect_tree_intersects_tree(n2->right_node, n1) )
		return LW_TRUE;
		else
		return LW_FALSE;
		}
		else
		{
		if( rect_tree_intersects_tree(n1->left_node, n2) \|\| rect_tree_intersects_tree(n1->right_node, n2) )
		return LW_TRUE;
		else
		return LW_FALSE;
		}
		}
		}
		else
		{
		return LW_FALSE;
		}
		}


		/**
		* Create a new leaf node, calculating a measure value for each point on the
		* edge and storing pointers back to the end points for later.
		*/
		RECT_NODE* rect_node_leaf_new(const POINTARRAY *pa, int i)
		{
		POINT2D p1, p2;
		RECT_NODE *node;

		p1 = (POINT2D*)getPoint_internal(pa, i);
		p2 = (POINT2D*)getPoint_internal(pa, i+1);

		/* Zero length edge, doesn't get a node */
		if( FP_EQUALS(p1->x, p2->x) && FP_EQUALS(p1->y, p2->y) )
		return NULL;

		node = lwalloc(sizeof(RECT_NODE));
		node->p1 = p1;
		node->p2 = p2;
		node->xmin = FP_MIN(p1->x,p2->x);
		node->xmax = FP_MAX(p1->x,p2->x);
		node->ymin = FP_MIN(p1->y,p2->y);
		node->ymax = FP_MAX(p1->y,p2->y);
		node->left_node = NULL;
		node->right_node = NULL;
		return node;
		}

		/**
		* Create a new internal node, calculating the new measure range for the node,
		* and storing pointers to the child nodes.
		*/
		RECT_NODE* rect_node_internal_new(RECT_NODE left_node, RECT_NODE right_node)
		{
		RECT_NODE *node = lwalloc(sizeof(RECT_NODE));
		node->p1 = NULL;
		node->p2 = NULL;
		node->xmin = FP_MIN(left_node->xmin, right_node->xmin);
		node->xmax = FP_MAX(left_node->xmax, right_node->xmax);
		node->ymin = FP_MIN(left_node->ymin, right_node->ymin);
		node->ymax = FP_MAX(left_node->ymax, right_node->ymax);
		node->left_node = left_node;
		node->right_node = right_node;
		return node;
		}

		/**
		* Build a tree of nodes from a point array, one node per edge, and each
		* with an associated measure range along a one-dimensional space. We
		* can then search that space as a range tree.
		*/
		RECT_NODE* rect_tree_new(const POINTARRAY *pa)
		{
		int num_edges, num_children, num_parents;
		int i, j;
		RECT_NODE **nodes;
		RECT_NODE *node;
		RECT_NODE *tree;

		if( pa->npoints < 2 )
		{
		return NULL;
		}

		/*
		** First create a flat list of nodes, one per edge.
		** For each vertex, transform into our one-dimensional measure.
		** Hopefully, when projected, the points turn into a fairly
		** uniformly distributed collection of measures.
		*/
		num_edges = pa->npoints - 1;
		nodes = lwalloc(sizeof(RECT_NODE) pa->npoints);
		j = 0;
		for( i = 0; i < num_edges; i++ )
		{
		node = rect_node_leaf_new(pa, i);
		if( node ) /* Not zero length? */
		{
		nodes[j] = node;
		j++;
		}
		}

		/*
		** If we sort the nodelist first, we'll get a more balanced tree
		** in the end, but at the cost of sorting. For now, we just
		** build the tree knowing that point arrays tend to have a
		** reasonable amount of sorting already.
		*/

		num_children = j;
		num_parents = num_children / 2;
		while( num_parents > 0 )
		{
		j = 0;
		while( j < num_parents )
		{
		/*
		** Each new parent includes pointers to the children, so even though
		** we are over-writing their place in the list, we still have references
		** to them via the tree.
		*/
		nodes[j] = rect_node_internal_new(nodes[2j], nodes[(2j)+1]);
		j++;
		}
		/* Odd number of children, just copy the last node up a level */
		if( num_children % 2 )
		{
		nodes[j] = nodes[num_children - 1];
		num_parents++;
		}
		num_children = num_parents;
		num_parents = num_children / 2;
		}

		/* Take a reference to the head of the tree*/
		tree = nodes[0];

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

		return tree;

		}