g_serialized.o \
g_util.o \
lwgeodetic.o \
- lwspheroid.o
+ lwspheroid.o \
+ lwtree.o
SA_HEADERS = \
liblwgeom.h \
** 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;
** 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;
/* 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;
}
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;
}
/* 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;
/* This should never happen! */
return SEG_ERROR;
-
}
/**
/* 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);
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 {
GEOGRAPHIC_POINT end;
} GEOGRAPHIC_EDGE;
+/**
+* Holder for sorting points in distance algorithm
+*/
+typedef struct
+{
+ double measure;
+ uint32 index;
+} DISTANCE_ORDER;
+
/**
* Conversion functions
*/
--- /dev/null
+#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[2*j], nodes[(2*j)+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;
+
+}
+
--- /dev/null
+#include "libgeom.h"
+
+
+/**
+* Note that p1 and p2 are pointers into an independent POINTARRAY, do not free them.
+*/
+typedef struct rect_node
+{
+ double xmin;
+ double xmax;
+ double ymin;
+ double ymax;
+ struct rect_node *left_node;
+ struct rect_node *right_node;
+ POINT2D *p1;
+ POINT2D *p2;
+} RECT_NODE;
+
+int rect_tree_contains_point(const RECT_NODE *tree, const POINT2D *pt, int *on_boundary);
+int rect_tree_intersects_tree(const RECT_NODE *tree1, const RECT_NODE *tree2);
+void rect_tree_free(RECT_NODE *tree);
+RECT_NODE* rect_node_leaf_new(const POINTARRAY *pa, int i);
+RECT_NODE* rect_node_internal_new(RECT_NODE *left_node, RECT_NODE *right_node);
+RECT_NODE* rect_tree_new(const POINTARRAY *pa);
thevalue = theP.x-(k*theP.y);
list1[t].themeasure=thevalue;
list1[t].pnr=t;
- /*lwnotice("1pnr=%d, pnr=%d, themeasure=%f, themeasure=%f",t,list1[t].pnr,thevalue,list1[t].themeasure );*/
+ //lwnotice("l1 %d, measure=%f",t,thevalue );
}
for (t=0; t<n2; t++) /*for each segment in L2*/
{
thevalue = theP.x-(k*theP.y);
list2[t].themeasure=thevalue;
list2[t].pnr=t;
+ //lwnotice("l2 %d, measure=%f",t,thevalue );
}
c1m = c1.x-(k*c1.y);
c2m = c2.x-(k*c2.y);
projects / postgis / commitdiff