/src/postgis/liblwgeom/gbox.c

Source
/**********************************************************************
 *
 * 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 <float.h>
#include <stdlib.h>
#include <math.h>

typedef struct {
  double x;
  double y;
  double z;
  double m;
  double w;
} NURBS_BBOX_HPOINT;

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_within_2d(const GBOX *g1, const GBOX *g2)
{
  if ( ( g1->xmin < g2->xmin ) || ( g1->xmax > g2->xmax ) ||
       ( g1->ymin < g2->ymin ) || ( g1->ymax > g2->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;
}

/**
* Warning, this function is only good for x/y/z boxes, used
* in unit testing of geodetic box generation.
*/
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 = {0};
  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;
}

/**
 * Projects a homogeneous control point by dividing its weighted coordinates
 * by w, failing when w is zero and the Cartesian point is undefined.
 */
static int
lwnurbscurve_hpoint_project(const NURBS_BBOX_HPOINT *hpoint, POINT4D *point)
{
  if (hpoint->w == 0.0)
    return LW_FAILURE;

  point->x = hpoint->x / hpoint->w;
  point->y = hpoint->y / hpoint->w;
  point->z = hpoint->z / hpoint->w;
  point->m = hpoint->m / hpoint->w;
  return LW_SUCCESS;
}

/**
 * Interpolates in homogeneous space, as required by rational de Casteljau
 * subdivision and Boehm knot insertion; interpolating projected points would
 * describe a different curve.
 */
static NURBS_BBOX_HPOINT
lwnurbscurve_hpoint_lerp(const NURBS_BBOX_HPOINT *a, const NURBS_BBOX_HPOINT *b, double t)
{
  NURBS_BBOX_HPOINT out;
  double s = 1.0 - t;

  out.x = s * a->x + t * b->x;
  out.y = s * a->y + t * b->y;
  out.z = s * a->z + t * b->z;
  out.m = s * a->m + t * b->m;
  out.w = s * a->w + t * b->w;
  return out;
}

/**
 * Loads one control point in homogeneous coordinates, applying the implicit
 * weight 1.0 for polynomial curves so Bezier extraction can use the same
 * arithmetic for rational and non-rational spans.
 */
static void
lwnurbscurve_get_hpoint(const LWNURBSCURVE *curve, uint32_t index, NURBS_BBOX_HPOINT *hpoint)
{
  POINT4D point;
  double weight = (curve->weights && index < curve->nweights) ? curve->weights[index] : 1.0;

  getPoint4d_p(curve->points, index, &point);
  hpoint->x = point.x * weight;
  hpoint->y = point.y * weight;
  hpoint->z = FLAGS_GET_Z(curve->flags) ? point.z * weight : 0.0;
  hpoint->m = FLAGS_GET_M(curve->flags) ? point.m * weight : 0.0;
  hpoint->w = weight;
}

/**
 * Expands a Cartesian box in each enabled dimension, keeping dimensionality
 * decisions out of the projection and subdivision helpers.
 */
static void
lwnurbscurve_add_point_to_gbox(const POINT4D *point, GBOX *gbox)
{
  gbox->xmin = FP_MIN(gbox->xmin, point->x);
  gbox->xmax = FP_MAX(gbox->xmax, point->x);
  gbox->ymin = FP_MIN(gbox->ymin, point->y);
  gbox->ymax = FP_MAX(gbox->ymax, point->y);

  if (FLAGS_GET_Z(gbox->flags))
  {
    gbox->zmin = FP_MIN(gbox->zmin, point->z);
    gbox->zmax = FP_MAX(gbox->zmax, point->z);
  }

  if (FLAGS_GET_M(gbox->flags))
  {
    gbox->mmin = FP_MIN(gbox->mmin, point->m);
    gbox->mmax = FP_MAX(gbox->mmax, point->m);
  }
}

/**
 * Projects and merges one homogeneous control point, centralizing zero-weight
 * rejection before any coordinate reaches a box accumulator.
 */
static int
lwnurbscurve_add_hpoint_to_gbox(const NURBS_BBOX_HPOINT *hpoint, GBOX *gbox)
{
  POINT4D point;

  if (lwnurbscurve_hpoint_project(hpoint, &point) == LW_FAILURE)
    return LW_FAILURE;

  lwnurbscurve_add_point_to_gbox(&point, gbox);
  return LW_SUCCESS;
}

/**
 * Initializes an empty double-precision accumulator for the requested
 * dimensions, using inverted extrema for later additions. Float rounding is
 * deferred until the subdivision stop test.
 */
static void
lwnurbscurve_gbox_init(GBOX *gbox, lwflags_t flags)
{
  gbox->flags = flags;
  gbox->xmin = gbox->ymin = gbox->zmin = gbox->mmin = DBL_MAX;
  gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = -DBL_MAX;
}

/**
 * Bounds the projected control-point hull of a positive-weight rational Bezier
 * span, providing the conservative outer bound used during subdivision.
 */
static int
lwnurbscurve_bezier_hull_gbox(const NURBS_BBOX_HPOINT *points, uint32_t npoints, lwflags_t flags, GBOX *gbox)
{
  uint32_t i;

  lwnurbscurve_gbox_init(gbox, flags);
  for (i = 0; i < npoints; i++)
  {
    if (lwnurbscurve_add_hpoint_to_gbox(&points[i], gbox) == LW_FAILURE)
      return LW_FAILURE;
  }

  return LW_SUCCESS;
}

/**
 * Bounds the projected endpoints of a Bezier span, providing an inner
 * coordinate-range bound to compare with the control hull's outer bound.
 */
static int
lwnurbscurve_bezier_endpoint_gbox(const NURBS_BBOX_HPOINT *points, uint32_t degree, lwflags_t flags, GBOX *gbox)
{
  lwnurbscurve_gbox_init(gbox, flags);

  if (lwnurbscurve_add_hpoint_to_gbox(&points[0], gbox) == LW_FAILURE ||
      lwnurbscurve_add_hpoint_to_gbox(&points[degree], gbox) == LW_FAILURE)
    return LW_FAILURE;

  return LW_SUCCESS;
}

/**
 * Tests whether two boxes round outward to the same float bounds, so refinement
 * stops once it can no longer change the serialized bounding box. The caller
 * supplies boxes with identical dimensional flags.
 */
static int
lwnurbscurve_gbox_same_float(const GBOX *a, const GBOX *b)
{
  if (next_float_down(a->xmin) != next_float_down(b->xmin) || next_float_up(a->xmax) != next_float_up(b->xmax) ||
      next_float_down(a->ymin) != next_float_down(b->ymin) || next_float_up(a->ymax) != next_float_up(b->ymax))
    return LW_FALSE;

  if (FLAGS_GET_Z(a->flags) &&
      (next_float_down(a->zmin) != next_float_down(b->zmin) || next_float_up(a->zmax) != next_float_up(b->zmax)))
    return LW_FALSE;

  if (FLAGS_GET_M(a->flags) &&
      (next_float_down(a->mmin) != next_float_down(b->mmin) || next_float_up(a->mmax) != next_float_up(b->mmax)))
    return LW_FALSE;

  return LW_TRUE;
}

/**
 * Detects when floating-point subdivision no longer changes a control net.
 * Exact equality is intentional: any changed component permits more progress;
 * otherwise recursion stops with the current conservative hull.
 */
static int
lwnurbscurve_hpoints_same(const NURBS_BBOX_HPOINT *a, const NURBS_BBOX_HPOINT *b, uint32_t npoints)
{
  uint32_t i;

  for (i = 0; i < npoints; i++)
  {
    if (a[i].x != b[i].x || a[i].y != b[i].y || a[i].z != b[i].z || a[i].m != b[i].m || a[i].w != b[i].w)
      return LW_FALSE;
  }

  return LW_TRUE;
}

/**
 * Splits a homogeneous Bezier control net at t = 1/2, producing two control
 * nets that represent the same curve on their respective half-intervals. For
 * positive weights, each child's projected control hull is contained in the
 * parent's hull, so subdivision can only tighten the conservative outer box.
 */
static void
lwnurbscurve_bezier_split(const NURBS_BBOX_HPOINT *points,
        uint32_t degree,
        NURBS_BBOX_HPOINT *left,
        NURBS_BBOX_HPOINT *right)
{
  uint32_t i, r;
  NURBS_BBOX_HPOINT tmp[degree + 1];

  for (i = 0; i <= degree; i++)
    tmp[i] = points[i];

  left[0] = tmp[0];
  right[degree] = tmp[degree];

  for (r = 1; r <= degree; r++)
  {
    for (i = 0; i <= degree - r; i++)
      tmp[i] = lwnurbscurve_hpoint_lerp(&tmp[i], &tmp[i + 1], 0.5);

    left[r] = tmp[0];
    right[degree - r] = tmp[degree - r];
  }
}

/**
 * Recursively bounds one rational Bezier span. Its projected control hull is a
 * conservative outer box, while its endpoints give an inner box. Subdivision
 * narrows the two child hulls and adds their shared split point to the endpoint
 * sample set. Once inner and outer boxes outward-round to the same float box,
 * the stored GSERIALIZED box is fixed.
 *
 * This helper owns the stop test, splitting, and recursion. In particular, it
 * does not merge an oversized parent hull before descending: only terminal
 * child hulls are merged into gbox.
 */
static int
lwnurbscurve_add_bezier_span_gbox(const NURBS_BBOX_HPOINT *points,
          uint32_t degree,
          lwflags_t flags,
          uint32_t depth,
          GBOX *gbox)
{
  GBOX hull_gbox, endpoint_gbox;
  NURBS_BBOX_HPOINT left[degree + 1];
  NURBS_BBOX_HPOINT right[degree + 1];

  if (lwnurbscurve_bezier_hull_gbox(points, degree + 1, flags, &hull_gbox) == LW_FAILURE ||
      lwnurbscurve_bezier_endpoint_gbox(points, degree, flags, &endpoint_gbox) == LW_FAILURE)
    return LW_FAILURE;

  /* A non-finite control point coordinate or weight (NaN or +/-Inf, accepted
   * by WKT input the same way LINESTRING accepts them) makes both stop tests
   * below unusable: NaN never compares equal to itself, so neither
   * lwnurbscurve_gbox_same_float() nor lwnurbscurve_hpoints_same() can ever
   * report convergence, and subdivision would run to the full DBL_MANT_DIG
   * recursion depth on every remaining span. hull_gbox already bounds every
   * control point in this span, so stop here and merge it directly, the same
   * way ptarray_calculate_gbox_cartesian() lets a non-finite point propagate
   * straight into the box for LINESTRING and the other simple types instead
   * of iterating further. */
  if (!gbox_is_valid(&hull_gbox))
  {
    gbox_merge(&hull_gbox, gbox);
    return LW_SUCCESS;
  }

  /* Positive NURBS weights put each rational Bezier span inside the convex
   * hull of its projected control points. Subdivision tightens that hull
   * until further refinement cannot change the float bbox that PostGIS stores
   * in GSERIALIZED. If arithmetic stops changing the control net first, merge
   * the current hull; it is still conservative. */
  if (depth >= DBL_MANT_DIG || lwnurbscurve_gbox_same_float(&hull_gbox, &endpoint_gbox))
  {
    gbox_merge(&hull_gbox, gbox);
    return LW_SUCCESS;
  }

  /* Keep splitting behind this helper so every child is subjected to the same
   * conservative stop and fallback rules before anything is merged. */
  lwnurbscurve_bezier_split(points, degree, left, right);
  if (lwnurbscurve_hpoints_same(left, points, degree + 1) || lwnurbscurve_hpoints_same(right, points, degree + 1))
  {
    gbox_merge(&hull_gbox, gbox);
    return LW_SUCCESS;
  }

  if (lwnurbscurve_add_bezier_span_gbox(left, degree, flags, depth + 1, gbox) == LW_FAILURE ||
      lwnurbscurve_add_bezier_span_gbox(right, degree, flags, depth + 1, gbox) == LW_FAILURE)
    return LW_FAILURE;

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

/**
 * Computes a conservative Cartesian box for the whole NURBS curve by
 * subdividing Bezier spans, which can tighten the box beyond what the original
 * control points or endpoints show. It submits each complete span to
 * lwnurbscurve_add_bezier_span_gbox(), which decides whether to split and owns
 * the recursive merge; splitting here would duplicate that policy.
 */
static int
lwnurbscurve_calculate_gbox_cartesian(const LWNURBSCURVE *curve, GBOX *gbox)
{
  uint32_t a, b, degree, i, j, m, multiplicity, npoints, r, s;
  int found_span = LW_FALSE;
  double domain_max, domain_min;
  double *knots;
  uint32_t nknots;

  if (!curve || !curve->points || curve->points->npoints == 0)
    return LW_FAILURE;

  degree = curve->degree;
  npoints = curve->points->npoints;
  knots = lwnurbscurve_get_or_generate_knots(curve, &nknots);
  if (!knots || nknots == 0)
    return LW_FAILURE;

  domain_min = knots[degree];
  domain_max = knots[npoints];
  if (domain_max <= domain_min || nknots != npoints + degree + 1)
  {
    lwfree(knots);
    return LW_FAILURE;
  }

  lwnurbscurve_gbox_init(gbox, lwflags(FLAGS_GET_Z(curve->flags), FLAGS_GET_M(curve->flags), 0));

  /* Stream the standard B-spline to Bezier decomposition with one current
   * and one next control net. This is Boehm knot insertion arranged so the
   * insertions for one knot block are applied only to the local control net,
   * preserving the exact span geometry while avoiding whole-curve knot
   * refinement and its quadratic memory churn. */
  {
    NURBS_BBOX_HPOINT bezier[degree + 1];
    NURBS_BBOX_HPOINT next_bezier[degree + 1];
    double alphas[degree ? degree : 1];

    for (i = 0; i <= degree; i++)
      lwnurbscurve_get_hpoint(curve, i, &bezier[i]);

    m = npoints + degree;
    a = degree;
    b = degree + 1;

    while (b < m)
    {
      uint32_t knot_block_start = b;

      /* Exact comparison is intentional: these values come directly from
       * the validated knot vector or its generated uniform replacement. */
      while (b < m && knots[b + 1] == knots[b])
        b++;

      multiplicity = b - knot_block_start + 1;
      if (multiplicity < degree)
      {
        double numerator = knots[b] - knots[a];

        /* Precompute the knot-insertion blend factors for this block.
         * A zero denominator would mean the validated active span has
         * collapsed, so fail instead of feeding NaNs to the bbox code. */
        for (j = degree; j > multiplicity; j--)
        {
          double denom = knots[a + j] - knots[a];

          if (denom == 0.0)
          {
            lwfree(knots);
            return LW_FAILURE;
          }

          alphas[j - multiplicity - 1] = numerator / denom;
        }

        r = degree - multiplicity;
        for (j = 1; j <= r; j++)
        {
          uint32_t save = r - j;

          /* Each pass raises the current knot multiplicity by one.
           * The rightmost point saved on every pass becomes the left
           * side of the next span's control net. */
          s = multiplicity + j;
          for (i = degree; i >= s; i--)
          {
            double alpha = alphas[i - s];
            bezier[i] = lwnurbscurve_hpoint_lerp(&bezier[i - 1], &bezier[i], alpha);
          }

          if (b < m)
            next_bezier[save] = bezier[degree];
        }
      }

      /* Only bbox spans inside the NURBS active parameter domain
       * [U[p], U[n]]. Unclamped explicit knot vectors can contain valid
       * knots before or after that domain, and those must not contribute. */
      if (knots[b] > knots[a] && knots[a] >= domain_min && knots[b] <= domain_max)
      {
        if (lwnurbscurve_add_bezier_span_gbox(bezier, degree, gbox->flags, 0, gbox) ==
            LW_FAILURE)
        {
          lwfree(knots);
          return LW_FAILURE;
        }
        found_span = LW_TRUE;
      }

      if (knots[b] >= domain_max || b >= m)
        break;

      /* Complete the next span from the saved right edge and the untouched
       * original control points that are entering the moving local window. */
      for (i = degree - FP_MIN(multiplicity, degree); i <= degree; i++)
        lwnurbscurve_get_hpoint(curve, b - degree + i, &next_bezier[i]);

      for (i = 0; i <= degree; i++)
        bezier[i] = next_bezier[i];

      a = b;
      b++;
    }
  }

  lwfree(knots);

  return found_span ? LW_SUCCESS : LW_FAILURE;
}

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 = {0};
  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;
}

/**
 * Calculate a Cartesian bounding box for an LWGEOM and store it in a GBOX.
 *
 * Computes the 2D/3D/M bounding box for the given geometry and writes the result
 * into the provided gbox. The function dispatches to type-specific helpers for
 * supported geometry types.
 *
 * @param lwgeom Geometry to compute the bounding box for. If NULL, the function returns LW_FAILURE.
 * @param gbox Destination GBOX to receive the computed bounding box (must be non-NULL).
 * @return LW_SUCCESS on successful computation and population of gbox; LW_FAILURE if lwgeom is NULL,
 *         if the geometry type is unsupported, or if a bounding box cannot be computed (e.g. empty NURBS).
 */
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);
  case NURBSCURVETYPE:
    return lwnurbscurve_calculate_gbox_cartesian((const LWNURBSCURVE *)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);
  }
}

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

Coverage Report

Created: 2026-08-31 06:50