gbox.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 2009 Paul Ramsey <pramsey@cleverelephant.ca>
 *
 **********************************************************************/
 
#include "liblwgeom_internal.h"
#include "lwgeodetic.h"
#include "lwgeom_log.h"
#include <stdlib.h>
#include <math.h>
 
 
GBOX* gbox_new(lwflags_t flags)
{
        GBOX *g = (GBOX*)lwalloc(sizeof(GBOX));
        gbox_init(g);
        g->flags = flags;
        return g;
}
 
void gbox_init(GBOX *gbox)
{
        memset(gbox, 0, sizeof(GBOX));
}
 
GBOX* gbox_clone(const GBOX *gbox)
{
        GBOX *g = lwalloc(sizeof(GBOX));
        memcpy(g, gbox, sizeof(GBOX));
        return g;
}
 
/* TODO to be removed */
BOX3D* box3d_from_gbox(const GBOX *gbox)
{
        BOX3D *b;
        assert(gbox);
 
        b = lwalloc(sizeof(BOX3D));
 
        b->xmin = gbox->xmin;
        b->xmax = gbox->xmax;
        b->ymin = gbox->ymin;
        b->ymax = gbox->ymax;
 
        if ( FLAGS_GET_Z(gbox->flags) )
        {
                b->zmin = gbox->zmin;
                b->zmax = gbox->zmax;
        }
        else
        {
                b->zmin = b->zmax = 0.0;
        }
 
        b->srid = SRID_UNKNOWN;
        return b;
}
 
/* TODO to be removed */
GBOX* box3d_to_gbox(const BOX3D *b3d)
{
        GBOX *b;
        assert(b3d);
 
        b = lwalloc(sizeof(GBOX));
 
        b->xmin = b3d->xmin;
        b->xmax = b3d->xmax;
        b->ymin = b3d->ymin;
        b->ymax = b3d->ymax;
        b->zmin = b3d->zmin;
        b->zmax = b3d->zmax;
 
        return b;
}
 
void gbox_expand(GBOX *g, double d)
{
        g->xmin -= d;
        g->xmax += d;
        g->ymin -= d;
        g->ymax += d;
        if (FLAGS_GET_Z(g->flags) || FLAGS_GET_GEODETIC(g->flags))
        {
                g->zmin -= d;
                g->zmax += d;
        }
        if (FLAGS_GET_M(g->flags))
        {
                g->mmin -= d;
                g->mmax += d;
        }
}
 
void gbox_expand_xyzm(GBOX *g, double dx, double dy, double dz, double dm)
{
        g->xmin -= dx;
        g->xmax += dx;
        g->ymin -= dy;
        g->ymax += dy;
 
        if (FLAGS_GET_Z(g->flags))
        {
                g->zmin -= dz;
                g->zmax += dz;
        }
 
        if (FLAGS_GET_M(g->flags))
        {
                g->mmin -= dm;
                g->mmax += dm;
        }
}
 
int gbox_union(const GBOX *g1, const GBOX *g2, GBOX *gout)
{
        if ( ( ! g1 ) && ( ! g2 ) )
                return LW_FALSE;
        else if (!g1)
        {
                memcpy(gout, g2, sizeof(GBOX));
                return LW_TRUE;
        }
        else if (!g2)
        {
                memcpy(gout, g1, sizeof(GBOX));
                return LW_TRUE;
        }
 
        gout->flags = g1->flags;
 
        gout->xmin = FP_MIN(g1->xmin, g2->xmin);
        gout->xmax = FP_MAX(g1->xmax, g2->xmax);
 
        gout->ymin = FP_MIN(g1->ymin, g2->ymin);
        gout->ymax = FP_MAX(g1->ymax, g2->ymax);
 
        gout->zmin = FP_MIN(g1->zmin, g2->zmin);
        gout->zmax = FP_MAX(g1->zmax, g2->zmax);
 
        return LW_TRUE;
}
 
int gbox_same(const GBOX *g1, const GBOX *g2)
{
        if (FLAGS_GET_ZM(g1->flags) != FLAGS_GET_ZM(g2->flags))
                return LW_FALSE;
 
        if (!gbox_same_2d(g1, g2)) return LW_FALSE;
 
        if (FLAGS_GET_Z(g1->flags) && (g1->zmin != g2->zmin || g1->zmax != g2->zmax))
                return LW_FALSE;
        if (FLAGS_GET_M(g1->flags) && (g1->mmin != g2->mmin || g1->mmax != g2->mmax))
                return LW_FALSE;
 
        return LW_TRUE;
}
 
int gbox_same_2d(const GBOX *g1, const GBOX *g2)
{
    if (g1->xmin == g2->xmin && g1->ymin == g2->ymin &&
        g1->xmax == g2->xmax && g1->ymax == g2->ymax)
                return LW_TRUE;
        return LW_FALSE;
}
 
int gbox_same_2d_float(const GBOX *g1, const GBOX *g2)
{
  if  ((g1->xmax == g2->xmax || next_float_up(g1->xmax)   == next_float_up(g2->xmax))   &&
       (g1->ymax == g2->ymax || next_float_up(g1->ymax)   == next_float_up(g2->ymax))   &&
       (g1->xmin == g2->xmin || next_float_down(g1->xmin) == next_float_down(g1->xmin)) &&
       (g1->ymin == g2->ymin || next_float_down(g2->ymin) == next_float_down(g2->ymin)))
      return LW_TRUE;
  return LW_FALSE;
}
 
int gbox_is_valid(const GBOX *gbox)
{
        /* X */
        if ( ! isfinite(gbox->xmin) || isnan(gbox->xmin) ||
             ! isfinite(gbox->xmax) || isnan(gbox->xmax) )
                return LW_FALSE;
 
        /* Y */
        if ( ! isfinite(gbox->ymin) || isnan(gbox->ymin) ||
             ! isfinite(gbox->ymax) || isnan(gbox->ymax) )
                return LW_FALSE;
 
        /* Z */
        if ( FLAGS_GET_GEODETIC(gbox->flags) || FLAGS_GET_Z(gbox->flags) )
        {
                if ( ! isfinite(gbox->zmin) || isnan(gbox->zmin) ||
                     ! isfinite(gbox->zmax) || isnan(gbox->zmax) )
                        return LW_FALSE;
        }
 
        /* M */
        if ( FLAGS_GET_M(gbox->flags) )
        {
                if ( ! isfinite(gbox->mmin) || isnan(gbox->mmin) ||
                     ! isfinite(gbox->mmax) || isnan(gbox->mmax) )
                        return LW_FALSE;
        }
 
        return LW_TRUE;
}
 
int gbox_merge_point3d(const POINT3D *p, GBOX *gbox)
{
        if ( gbox->xmin > p->x ) gbox->xmin = p->x;
        if ( gbox->ymin > p->y ) gbox->ymin = p->y;
        if ( gbox->zmin > p->z ) gbox->zmin = p->z;
        if ( gbox->xmax < p->x ) gbox->xmax = p->x;
        if ( gbox->ymax < p->y ) gbox->ymax = p->y;
        if ( gbox->zmax < p->z ) gbox->zmax = p->z;
        return LW_SUCCESS;
}
 
int gbox_init_point3d(const POINT3D *p, GBOX *gbox)
{
        gbox->xmin = gbox->xmax = p->x;
        gbox->ymin = gbox->ymax = p->y;
        gbox->zmin = gbox->zmax = p->z;
        return LW_SUCCESS;
}
 
int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
{
        if ( gbox->xmin > pt->x || gbox->ymin > pt->y || gbox->zmin > pt->z ||
             gbox->xmax < pt->x || gbox->ymax < pt->y || gbox->zmax < pt->z )
        {
                return LW_FALSE;
        }
        return LW_TRUE;
}
 
int gbox_merge(const GBOX *new_box, GBOX *merge_box)
{
        assert(merge_box);
 
        if ( FLAGS_GET_ZM(merge_box->flags) != FLAGS_GET_ZM(new_box->flags) )
                return LW_FAILURE;
 
        if ( new_box->xmin < merge_box->xmin) merge_box->xmin = new_box->xmin;
        if ( new_box->ymin < merge_box->ymin) merge_box->ymin = new_box->ymin;
        if ( new_box->xmax > merge_box->xmax) merge_box->xmax = new_box->xmax;
        if ( new_box->ymax > merge_box->ymax) merge_box->ymax = new_box->ymax;
 
        if ( FLAGS_GET_Z(merge_box->flags) || FLAGS_GET_GEODETIC(merge_box->flags) )
        {
                if ( new_box->zmin < merge_box->zmin) merge_box->zmin = new_box->zmin;
                if ( new_box->zmax > merge_box->zmax) merge_box->zmax = new_box->zmax;
        }
        if ( FLAGS_GET_M(merge_box->flags) )
        {
                if ( new_box->mmin < merge_box->mmin) merge_box->mmin = new_box->mmin;
                if ( new_box->mmax > merge_box->mmax) merge_box->mmax = new_box->mmax;
        }
 
        return LW_SUCCESS;
}
 
int gbox_overlaps(const GBOX *g1, const GBOX *g2)
{
 
        /* Make sure our boxes are consistent */
        if ( FLAGS_GET_GEODETIC(g1->flags) != FLAGS_GET_GEODETIC(g2->flags) )
                lwerror("gbox_overlaps: cannot compare geodetic and non-geodetic boxes");
 
        /* Check X/Y first */
        if ( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
             g1->xmin > g2->xmax || g1->ymin > g2->ymax )
                return LW_FALSE;
 
        /* Deal with the geodetic case special: we only compare the geodetic boxes (x/y/z) */
        /* Never the M dimension */
        if ( FLAGS_GET_GEODETIC(g1->flags) && FLAGS_GET_GEODETIC(g2->flags) )
        {
                if ( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
                        return LW_FALSE;
                else
                        return LW_TRUE;
        }
 
        /* If both geodetic or both have Z, check Z */
        if ( FLAGS_GET_Z(g1->flags) && FLAGS_GET_Z(g2->flags) )
        {
                if ( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
                        return LW_FALSE;
        }
 
        /* If both have M, check M */
        if ( FLAGS_GET_M(g1->flags) && FLAGS_GET_M(g2->flags) )
        {
                if ( g1->mmax < g2->mmin || g1->mmin > g2->mmax )
                        return LW_FALSE;
        }
 
        return LW_TRUE;
}
 
int
gbox_overlaps_2d(const GBOX *g1, const GBOX *g2)
{
 
        /* Make sure our boxes are consistent */
        if ( FLAGS_GET_GEODETIC(g1->flags) != FLAGS_GET_GEODETIC(g2->flags) )
                lwerror("gbox_overlaps: cannot compare geodetic and non-geodetic boxes");
 
        /* Check X/Y first */
        if ( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
             g1->xmin > g2->xmax || g1->ymin > g2->ymax )
                return LW_FALSE;
 
        return LW_TRUE;
}
 
int
gbox_contains_2d(const GBOX *g1, const GBOX *g2)
{
        if ( ( g2->xmin < g1->xmin ) || ( g2->xmax > g1->xmax ) ||
             ( g2->ymin < g1->ymin ) || ( g2->ymax > g1->ymax ) )
        {
                return LW_FALSE;
        }
        return LW_TRUE;
}
 
int
gbox_contains_point2d(const GBOX *g, const POINT2D *p)
{
        if ( ( g->xmin <= p->x ) && ( g->xmax >= p->x ) &&
             ( g->ymin <= p->y ) && ( g->ymax >= p->y ) )
        {
                return LW_TRUE;
        }
        return LW_FALSE;
}
 
GBOX* gbox_from_string(const char *str)
{
        const char *ptr = str;
        char *nextptr;
        char *gbox_start = strstr(str, "GBOX((");
        GBOX *gbox = gbox_new(lwflags(0,0,1));
        if ( ! gbox_start ) return NULL; /* No header found */
        ptr += 6;
        gbox->xmin = strtod(ptr, &nextptr);
        if ( ptr == nextptr ) return NULL; /* No double found */
        ptr = nextptr + 1;
        gbox->ymin = strtod(ptr, &nextptr);
        if ( ptr == nextptr ) return NULL; /* No double found */
        ptr = nextptr + 1;
        gbox->zmin = strtod(ptr, &nextptr);
        if ( ptr == nextptr ) return NULL; /* No double found */
        ptr = nextptr + 3;
        gbox->xmax = strtod(ptr, &nextptr);
        if ( ptr == nextptr ) return NULL; /* No double found */
        ptr = nextptr + 1;
        gbox->ymax = strtod(ptr, &nextptr);
        if ( ptr == nextptr ) return NULL; /* No double found */
        ptr = nextptr + 1;
        gbox->zmax = strtod(ptr, &nextptr);
        if ( ptr == nextptr ) return NULL; /* No double found */
        return gbox;
}
 
char* gbox_to_string(const GBOX *gbox)
{
        const size_t sz = 138;
        char *str = NULL;
 
        if ( ! gbox )
                return lwstrdup("NULL POINTER");
 
        str = (char*)lwalloc(sz);
 
        if ( FLAGS_GET_GEODETIC(gbox->flags) )
        {
                snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->xmax, gbox->ymax, gbox->zmax);
                return str;
        }
        if ( FLAGS_GET_Z(gbox->flags) && FLAGS_GET_M(gbox->flags) )
        {
                snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->mmin, gbox->xmax, gbox->ymax, gbox->zmax, gbox->mmax);
                return str;
        }
        if ( FLAGS_GET_Z(gbox->flags) )
        {
                snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->xmax, gbox->ymax, gbox->zmax);
                return str;
        }
        if ( FLAGS_GET_M(gbox->flags) )
        {
                snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->mmin, gbox->xmax, gbox->ymax, gbox->mmax);
                return str;
        }
        snprintf(str, sz, "GBOX((%.8g,%.8g),(%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->xmax, gbox->ymax);
        return str;
}
 
GBOX* gbox_copy(const GBOX *box)
{
        GBOX *copy = (GBOX*)lwalloc(sizeof(GBOX));
        memcpy(copy, box, sizeof(GBOX));
        return copy;
}
 
void gbox_duplicate(const GBOX *original, GBOX *duplicate)
{
        assert(duplicate);
        assert(original);
        memcpy(duplicate, original, sizeof(GBOX));
}
 
size_t gbox_serialized_size(lwflags_t flags)
{
        if (FLAGS_GET_GEODETIC(flags))
                return 6 * sizeof(float);
        else
                return 2 * FLAGS_NDIMS(flags) * sizeof(float);
}
 
 
/* ********************************************************************************
** Compute cartesian bounding GBOX boxes from LWGEOM.
*/
 
int lw_arc_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox)
{
        POINT2D xmin, ymin, xmax, ymax;
        POINT2D C;
        int A2_side;
        double radius_A;
 
        LWDEBUG(2, "lw_arc_calculate_gbox_cartesian_2d called.");
 
        radius_A = lw_arc_center(A1, A2, A3, &C);
 
        /* Negative radius signals straight line, p1/p2/p3 are collinear */
        if (radius_A < 0.0)
        {
        gbox->xmin = FP_MIN(A1->x, A3->x);
        gbox->ymin = FP_MIN(A1->y, A3->y);
        gbox->xmax = FP_MAX(A1->x, A3->x);
        gbox->ymax = FP_MAX(A1->y, A3->y);
            return LW_SUCCESS;
        }
 
        /* Matched start/end points imply circle */
        if ( A1->x == A3->x && A1->y == A3->y )
        {
                gbox->xmin = C.x - radius_A;
                gbox->ymin = C.y - radius_A;
                gbox->xmax = C.x + radius_A;
                gbox->ymax = C.y + radius_A;
                return LW_SUCCESS;
        }
 
        /* First approximation, bounds of start/end points */
    gbox->xmin = FP_MIN(A1->x, A3->x);
    gbox->ymin = FP_MIN(A1->y, A3->y);
    gbox->xmax = FP_MAX(A1->x, A3->x);
    gbox->ymax = FP_MAX(A1->y, A3->y);
 
        /* Create points for the possible extrema */
        xmin.x = C.x - radius_A;
        xmin.y = C.y;
        ymin.x = C.x;
        ymin.y = C.y - radius_A;
        xmax.x = C.x + radius_A;
        xmax.y = C.y;
        ymax.x = C.x;
        ymax.y = C.y + radius_A;
 
        /* Divide the circle into two parts, one on each side of a line
           joining p1 and p3. The circle extrema on the same side of that line
           as p2 is on, are also the extrema of the bbox. */
 
        A2_side = lw_segment_side(A1, A3, A2);
 
        if ( A2_side == lw_segment_side(A1, A3, &xmin) )
                gbox->xmin = xmin.x;
 
        if ( A2_side == lw_segment_side(A1, A3, &ymin) )
                gbox->ymin = ymin.y;
 
        if ( A2_side == lw_segment_side(A1, A3, &xmax) )
                gbox->xmax = xmax.x;
 
        if ( A2_side == lw_segment_side(A1, A3, &ymax) )
                gbox->ymax = ymax.y;
 
        return LW_SUCCESS;
}
 
 
static int lw_arc_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox)
{
        int rv;
 
        LWDEBUG(2, "lw_arc_calculate_gbox_cartesian called.");
 
        rv = lw_arc_calculate_gbox_cartesian_2d((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, gbox);
    gbox->zmin = FP_MIN(p1->z, p3->z);
    gbox->mmin = FP_MIN(p1->m, p3->m);
    gbox->zmax = FP_MAX(p1->z, p3->z);
    gbox->mmax = FP_MAX(p1->m, p3->m);
        return rv;
}
 
static void
ptarray_calculate_gbox_cartesian_2d(const POINTARRAY *pa, GBOX *gbox)
{
        const POINT2D *p = getPoint2d_cp(pa, 0);
 
        gbox->xmax = gbox->xmin = p->x;
        gbox->ymax = gbox->ymin = p->y;
 
        for (uint32_t i = 1; i < pa->npoints; i++)
        {
                p = getPoint2d_cp(pa, i);
                gbox->xmin = FP_MIN(gbox->xmin, p->x);
                gbox->xmax = FP_MAX(gbox->xmax, p->x);
                gbox->ymin = FP_MIN(gbox->ymin, p->y);
                gbox->ymax = FP_MAX(gbox->ymax, p->y);
        }
}
 
/* Works with X/Y/Z. Needs to be adjusted after if X/Y/M was required */
static void
ptarray_calculate_gbox_cartesian_3d(const POINTARRAY *pa, GBOX *gbox)
{
        const POINT3D *p = getPoint3d_cp(pa, 0);
 
        gbox->xmax = gbox->xmin = p->x;
        gbox->ymax = gbox->ymin = p->y;
        gbox->zmax = gbox->zmin = p->z;
 
        for (uint32_t i = 1; i < pa->npoints; i++)
        {
                p = getPoint3d_cp(pa, i);
                gbox->xmin = FP_MIN(gbox->xmin, p->x);
                gbox->xmax = FP_MAX(gbox->xmax, p->x);
                gbox->ymin = FP_MIN(gbox->ymin, p->y);
                gbox->ymax = FP_MAX(gbox->ymax, p->y);
                gbox->zmin = FP_MIN(gbox->zmin, p->z);
                gbox->zmax = FP_MAX(gbox->zmax, p->z);
        }
}
 
static void
ptarray_calculate_gbox_cartesian_4d(const POINTARRAY *pa, GBOX *gbox)
{
        const POINT4D *p = getPoint4d_cp(pa, 0);
 
        gbox->xmax = gbox->xmin = p->x;
        gbox->ymax = gbox->ymin = p->y;
        gbox->zmax = gbox->zmin = p->z;
        gbox->mmax = gbox->mmin = p->m;
 
        for (uint32_t i = 1; i < pa->npoints; i++)
        {
                p = getPoint4d_cp(pa, i);
                gbox->xmin = FP_MIN(gbox->xmin, p->x);
                gbox->xmax = FP_MAX(gbox->xmax, p->x);
                gbox->ymin = FP_MIN(gbox->ymin, p->y);
                gbox->ymax = FP_MAX(gbox->ymax, p->y);
                gbox->zmin = FP_MIN(gbox->zmin, p->z);
                gbox->zmax = FP_MAX(gbox->zmax, p->z);
                gbox->mmin = FP_MIN(gbox->mmin, p->m);
                gbox->mmax = FP_MAX(gbox->mmax, p->m);
        }
}
 
int
ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
{
        if (!pa || pa->npoints == 0)
                return LW_FAILURE;
        if (!gbox)
                return LW_FAILURE;
 
        int has_z = FLAGS_GET_Z(pa->flags);
        int has_m = FLAGS_GET_M(pa->flags);
        gbox->flags = lwflags(has_z, has_m, 0);
        LWDEBUGF(4, "ptarray_calculate_gbox Z: %d M: %d", has_z, has_m);
        int coordinates = 2 + has_z + has_m;
 
        switch (coordinates)
        {
        case 2:
        {
                ptarray_calculate_gbox_cartesian_2d(pa, gbox);
                break;
        }
        case 3:
        {
                if (has_z)
                {
                        ptarray_calculate_gbox_cartesian_3d(pa, gbox);
                }
                else
                {
                        double zmin = gbox->zmin;
                        double zmax = gbox->zmax;
                        ptarray_calculate_gbox_cartesian_3d(pa, gbox);
                        gbox->mmin = gbox->zmin;
                        gbox->mmax = gbox->zmax;
                        gbox->zmin = zmin;
                        gbox->zmax = zmax;
                }
                break;
        }
        default:
        {
                ptarray_calculate_gbox_cartesian_4d(pa, gbox);
                break;
        }
        }
        return LW_SUCCESS;
}
 
static int lwcircstring_calculate_gbox_cartesian(LWCIRCSTRING *curve, GBOX *gbox)
{
        GBOX tmp;
        POINT4D p1, p2, p3;
        uint32_t i;
 
        if (!curve) return LW_FAILURE;
        if (curve->points->npoints < 3) return LW_FAILURE;
 
        tmp.flags =
            lwflags(FLAGS_GET_Z(curve->flags), FLAGS_GET_M(curve->flags), 0);
 
        /* Initialize */
        gbox->xmin = gbox->ymin = gbox->zmin = gbox->mmin = FLT_MAX;
        gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = -1*FLT_MAX;
 
        for ( i = 2; i < curve->points->npoints; i += 2 )
        {
                getPoint4d_p(curve->points, i-2, &p1);
                getPoint4d_p(curve->points, i-1, &p2);
                getPoint4d_p(curve->points, i, &p3);
 
                if (lw_arc_calculate_gbox_cartesian(&p1, &p2, &p3, &tmp) == LW_FAILURE)
                        continue;
 
                gbox_merge(&tmp, gbox);
        }
 
        return LW_SUCCESS;
}
 
static int lwpoint_calculate_gbox_cartesian(LWPOINT *point, GBOX *gbox)
{
        if ( ! point ) return LW_FAILURE;
        return ptarray_calculate_gbox_cartesian( point->point, gbox );
}
 
static int lwline_calculate_gbox_cartesian(LWLINE *line, GBOX *gbox)
{
        if ( ! line ) return LW_FAILURE;
        return ptarray_calculate_gbox_cartesian( line->points, gbox );
}
 
static int lwtriangle_calculate_gbox_cartesian(LWTRIANGLE *triangle, GBOX *gbox)
{
        if ( ! triangle ) return LW_FAILURE;
        return ptarray_calculate_gbox_cartesian( triangle->points, gbox );
}
 
static int lwpoly_calculate_gbox_cartesian(LWPOLY *poly, GBOX *gbox)
{
        if ( ! poly ) return LW_FAILURE;
        if ( poly->nrings == 0 ) return LW_FAILURE;
        /* Just need to check outer ring */
        return ptarray_calculate_gbox_cartesian( poly->rings[0], gbox );
}
 
static int lwcollection_calculate_gbox_cartesian(LWCOLLECTION *coll, GBOX *gbox)
{
        GBOX subbox;
        uint32_t i;
        int result = LW_FAILURE;
        int first = LW_TRUE;
        assert(coll);
        if ( (coll->ngeoms == 0) || !gbox)
                return LW_FAILURE;
 
        subbox.flags = coll->flags;
 
        for ( i = 0; i < coll->ngeoms; i++ )
        {
                if ( lwgeom_calculate_gbox_cartesian((LWGEOM*)(coll->geoms[i]), &subbox) == LW_SUCCESS )
                {
                        /* Keep a copy of the sub-bounding box for later
                        if ( coll->geoms[i]->bbox )
                                lwfree(coll->geoms[i]->bbox);
                        coll->geoms[i]->bbox = gbox_copy(&subbox); */
                        if ( first )
                        {
                                gbox_duplicate(&subbox, gbox);
                                first = LW_FALSE;
                        }
                        else
                        {
                                gbox_merge(&subbox, gbox);
                        }
                        result = LW_SUCCESS;
                }
        }
        return result;
}
 
int lwgeom_calculate_gbox_cartesian(const LWGEOM *lwgeom, GBOX *gbox)
{
        if ( ! lwgeom ) return LW_FAILURE;
        LWDEBUGF(4, "lwgeom_calculate_gbox got type (%d) - %s", lwgeom->type, lwtype_name(lwgeom->type));
 
        switch (lwgeom->type)
        {
        case POINTTYPE:
                return lwpoint_calculate_gbox_cartesian((LWPOINT *)lwgeom, gbox);
        case LINETYPE:
                return lwline_calculate_gbox_cartesian((LWLINE *)lwgeom, gbox);
        case CIRCSTRINGTYPE:
                return lwcircstring_calculate_gbox_cartesian((LWCIRCSTRING *)lwgeom, gbox);
        case POLYGONTYPE:
                return lwpoly_calculate_gbox_cartesian((LWPOLY *)lwgeom, gbox);
        case TRIANGLETYPE:
                return lwtriangle_calculate_gbox_cartesian((LWTRIANGLE *)lwgeom, gbox);
        case COMPOUNDTYPE:
        case CURVEPOLYTYPE:
        case MULTIPOINTTYPE:
        case MULTILINETYPE:
        case MULTICURVETYPE:
        case MULTIPOLYGONTYPE:
        case MULTISURFACETYPE:
        case POLYHEDRALSURFACETYPE:
        case TINTYPE:
        case COLLECTIONTYPE:
                return lwcollection_calculate_gbox_cartesian((LWCOLLECTION *)lwgeom, gbox);
        }
        /* Never get here, please. */
        lwerror("unsupported type (%d) - %s", lwgeom->type, lwtype_name(lwgeom->type));
        return LW_FAILURE;
}
 
void gbox_float_round(GBOX *gbox)
{
        gbox->xmin = next_float_down(gbox->xmin);
        gbox->xmax = next_float_up(gbox->xmax);
 
        gbox->ymin = next_float_down(gbox->ymin);
        gbox->ymax = next_float_up(gbox->ymax);
 
        if ( FLAGS_GET_M(gbox->flags) )
        {
                gbox->mmin = next_float_down(gbox->mmin);
                gbox->mmax = next_float_up(gbox->mmax);
        }
 
        if ( FLAGS_GET_Z(gbox->flags) )
        {
                gbox->zmin = next_float_down(gbox->zmin);
                gbox->zmax = next_float_up(gbox->zmax);
        }
}
 
inline static uint64_t
uint64_interleave_2(uint64_t x, uint64_t y)
{
        x = (x | (x << 16)) & 0x0000FFFF0000FFFFULL;
        x = (x | (x << 8)) & 0x00FF00FF00FF00FFULL;
        x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0FULL;
        x = (x | (x << 2)) & 0x3333333333333333ULL;
        x = (x | (x << 1)) & 0x5555555555555555ULL;
 
        y = (y | (y << 16)) & 0x0000FFFF0000FFFFULL;
        y = (y | (y << 8)) & 0x00FF00FF00FF00FFULL;
        y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0FULL;
        y = (y | (y << 2)) & 0x3333333333333333ULL;
        y = (y | (y << 1)) & 0x5555555555555555ULL;
 
        return x | (y << 1);
}
 
/* Based on https://github.com/rawrunprotected/hilbert_curves Public Domain code */
inline static uint64_t
uint32_hilbert(uint32_t px, uint32_t py)
{
        uint64_t x = px;
        uint64_t y = py;
 
        uint64_t A, B, C, D;
 
        // Initial prefix scan round, prime with x and y
        {
                uint64_t a = x ^ y;
                uint64_t b = 0xFFFFFFFFULL ^ a;
                uint64_t c = 0xFFFFFFFFULL ^ (x | y);
                uint64_t d = x & (y ^ 0xFFFFFFFFULL);
 
                A = a | (b >> 1);
                B = (a >> 1) ^ a;
                C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
                D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
        }
 
        {
                uint64_t a = A;
                uint64_t b = B;
                uint64_t c = C;
                uint64_t d = D;
 
                A = ((a & (a >> 2)) ^ (b & (b >> 2)));
                B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
                C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
                D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
        }
 
        {
                uint64_t a = A;
                uint64_t b = B;
                uint64_t c = C;
                uint64_t d = D;
 
                A = ((a & (a >> 4)) ^ (b & (b >> 4)));
                B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
                C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
                D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
        }
 
        {
                uint64_t a = A;
                uint64_t b = B;
                uint64_t c = C;
                uint64_t d = D;
 
                A = ((a & (a >> 8)) ^ (b & (b >> 8)));
                B = ((a & (b >> 8)) ^ (b & ((a ^ b) >> 8)));
                C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
                D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
        }
 
        {
                uint64_t a = A;
                uint64_t b = B;
                uint64_t c = C;
                uint64_t d = D;
 
                C ^= ((a & (c >> 16)) ^ (b & (d >> 16)));
                D ^= ((b & (c >> 16)) ^ ((a ^ b) & (d >> 16)));
        }
 
        // Undo transformation prefix scan
        uint64_t a = C ^ (C >> 1);
        uint64_t b = D ^ (D >> 1);
 
        // Recover index bits
        uint64_t i0 = x ^ y;
        uint64_t i1 = b | (0xFFFFFFFFULL ^ (i0 | a));
 
        return uint64_interleave_2(i0, i1);
}
 
uint64_t
gbox_get_sortable_hash(const GBOX *g, const int32_t srid)
{
        union floatuint {
                uint32_t u;
                float f;
        };
 
        union floatuint x, y;
 
        /*
        * Since in theory the bitwise representation of an IEEE
        * float is sortable (exponents come before mantissa, etc)
        * we just copy the bits directly into an int and then
        * interleave those ints.
        */
        if (FLAGS_GET_GEODETIC(g->flags))
        {
                GEOGRAPHIC_POINT gpt;
                POINT3D p;
                p.x = (g->xmax + g->xmin) / 2.0;
                p.y = (g->ymax + g->ymin) / 2.0;
                p.z = (g->zmax + g->zmin) / 2.0;
                normalize(&p);
                cart2geog(&p, &gpt);
                /* We know range for geography, so build the curve taking it into account */
                x.f = 1.5 + gpt.lon / 512.0;
                y.f = 1.5 + gpt.lat / 256.0;
        }
        else
        {
                x.f = (g->xmax + g->xmin) / 2;
                y.f = (g->ymax + g->ymin) / 2;
                /*
                 * Tweak for popular SRID values: push floating point values into 1..2 range,
                 * a region where exponent is constant and thus Hilbert curve
                 * doesn't have compression artifact when X or Y value is close to 0.
                 * If someone has out of bounds value it will still expose the arifact but not crash.
                 * TODO: reconsider when we will have machinery to properly get bounds by SRID.
                 */
                if (srid == 3857 || srid == 3395)
                {
                        x.f = 1.5 + x.f / 67108864.0;
                        y.f = 1.5 + y.f / 67108864.0;
                }
                else if (srid == 4326)
                {
                        x.f = 1.5 + x.f / 512.0;
                        y.f = 1.5 + y.f / 256.0;
                }
        }
 
        return uint32_hilbert(y.u, x.u);
}