lwtree.c

Go to the documentation of this file.

/**********************************************************************
 *
 * 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"


static int rect_node_is_leaf(const RECT_NODE *node)
{
        return (node->p1 != NULL);
}

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 )
                                *on_boundary = LW_TRUE;
                        return (side < 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)
{
        LWDEBUGF(4,"n1 (%.9g %.9g,%.9g %.9g) vs n2 (%.9g %.9g,%.9g %.9g)",n1->xmin,n1->ymin,n1->xmax,n1->ymax,n2->xmin,n2->ymin,n2->xmax,n2->ymax);
        /* There can only be an edge intersection if the rectangles overlap */
        if ( ! ( FP_GT(n1->xmin, n2->xmax) || FP_GT(n2->xmin, n1->xmax) || FP_GT(n1->ymin, n2->ymax) || FP_GT(n2->ymin, n1->ymax) ) )
        {
                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 */
                        if ( lw_segment_intersects(n1->p1, n1->p2, n2->p1, n2->p2) )
                                return LW_TRUE;
                        else
                                return LW_FALSE;
                }
                else
                {
                        LWDEBUG(4,"  internal node found, recursing");
                        /* 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
        {
                LWDEBUG(4," no interaction found");
                return LW_FALSE;
        }
}


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;
}

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;
}

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;

}