Speed up the calculation of cartesian bbox

Closes https://github.com/postgis/postgis/pull/475
Closes https://trac.osgeo.org/postgis/ticket/4503

git-svn-id: http://svn.osgeo.org/postgis/trunk@17808 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/NEWS b/NEWS
index ae50678c1cf2f40dfc9f613c97b4707b0c88fc7f..e22bd6775572c564c52eda9399656bc067f736ed 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -11,6 +11,7 @@ Additional features enabled if you are running Proj6+ and PostgreSQL 12
   - #4495, Fix ST_SnapToGrid output having an outdated bbox (Raúl Marín)
   - #4496, Make ST_Simplify(TRIANGLE) collapse if requested (Raúl Marín)
   - #4501, Allow postgis_tiger_geocoder to be installable by non-super users (Regina Obe)
+  - #4503, Speed up the calculation of cartesian bbox (Raúl Marín)
 
 PostGIS 3.0.0alpha4
 2019/08/10

diff --git a/liblwgeom/gbox.c b/liblwgeom/gbox.c
index 31ea45bdec3c84a92feb3c6538a387d26a992206..141f3c2647e7817d61ae6caf2c8a3109f89ccde2 100644 (file)
--- a/liblwgeom/gbox.c
+++ b/liblwgeom/gbox.c
@@ -533,46 +533,114 @@ static int lw_arc_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2,
        return rv;
 }
 
-int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox )
+static void
+ptarray_calculate_gbox_cartesian_2d(const POINTARRAY *pa, GBOX *gbox)
 {
-       uint32_t i;
-       POINT4D p;
-       int has_z, has_m;
+       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;
 
-       if ( ! pa ) return LW_FAILURE;
-       if ( ! gbox ) return LW_FAILURE;
-       if ( pa->npoints < 1 ) return LW_FAILURE;
+       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;
 
-       has_z = FLAGS_GET_Z(pa->flags);
-       has_m = FLAGS_GET_M(pa->flags);
+       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;
 
-       getPoint4d_p(pa, 0, &p);
-       gbox->xmin = gbox->xmax = p.x;
-       gbox->ymin = gbox->ymax = p.y;
-       if ( has_z )
-               gbox->zmin = gbox->zmax = p.z;
-       if ( has_m )
-               gbox->mmin = gbox->mmax = p.m;
-
-       for ( i = 1 ; i < pa->npoints; i++ )
-       {
-               getPoint4d_p(pa, i, &p);
-               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);
-               if ( has_z )
+       switch (coordinates)
+       {
+       case 2:
+       {
+               ptarray_calculate_gbox_cartesian_2d(pa, gbox);
+               break;
+       }
+       case 3:
+       {
+               if (has_z)
                {
-                       gbox->zmin = FP_MIN(gbox->zmin, p.z);
-                       gbox->zmax = FP_MAX(gbox->zmax, p.z);
+                       ptarray_calculate_gbox_cartesian_3d(pa, gbox);
                }
-               if ( has_m )
+               else
                {
-                       gbox->mmin = FP_MIN(gbox->mmin, p.m);
-                       gbox->mmax = FP_MAX(gbox->mmax, p.m);
+                       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;
 }

diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 75f9fd6c297e934567f1c82be018b95b17a9db43..e12efdc79a285576975ee7aac7c5df123836c941 100644 (file)
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -953,22 +953,6 @@ extern POINT2D getPoint2d(const POINTARRAY *pa, uint32_t n);
  */
 extern int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point);
 
-/**
-* Returns a POINT3DZ pointer into the POINTARRAY serialized_ptlist,
-* suitable for reading from. This is very high performance
-* and declared const because you aren't allowed to muck with the
-* values, only read them.
-*/
-extern const POINT3DZ* getPoint3dz_cp(const POINTARRAY *pa, uint32_t n);
-
-/**
-* Returns a POINT4D pointer into the POINTARRAY serialized_ptlist,
-* suitable for reading from. This is very high performance
-* and declared const because you aren't allowed to muck with the
-* values, only read them.
-*/
-extern const POINT4D* getPoint4d_cp(const POINTARRAY *pa, uint32_t n);
-
 /*
  * set point N to the given value
  * NOTE that the pointarray can be of any

diff --git a/liblwgeom/lwgeom_api.c b/liblwgeom/lwgeom_api.c
index dfbcf18632e4d207626078eeb440de5ad6bced13..90e7a99f847826ec4bdc7435d6146479f2e1452b 100644 (file)
--- a/liblwgeom/lwgeom_api.c
+++ b/liblwgeom/lwgeom_api.c
@@ -362,46 +362,6 @@ getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
        return 1;
 }
 
-const POINT3DZ*
-getPoint3dz_cp(const POINTARRAY *pa, uint32_t n)
-{
-       if ( ! pa ) return 0;
-
-       if ( ! FLAGS_GET_Z(pa->flags) )
-       {
-               lwerror("getPoint3dz_cp: no Z coordinates in point array");
-               return 0; /*error */
-       }
-
-       if ( n>=pa->npoints )
-       {
-               lwerror("getPoint3dz_cp: point offset out of range");
-               return 0; /*error */
-       }
-
-       return (const POINT3DZ*)getPoint_internal(pa, n);
-}
-
-const POINT4D*
-getPoint4d_cp(const POINTARRAY* pa, uint32_t n)
-{
-       if (!pa) return 0;
-
-       if (!(FLAGS_GET_Z(pa->flags) && FLAGS_GET_M(pa->flags)))
-       {
-               lwerror("getPoint4d_cp: no Z and M coordinates in point array");
-               return 0; /*error */
-       }
-
-       if (n >= pa->npoints)
-       {
-               lwerror("getPoint4d_cp: point offset out of range");
-               return 0; /*error */
-       }
-
-       return (const POINT4D*)getPoint_internal(pa, n);
-}
-
 /*
  * set point N to the given value
  * NOTE that the pointarray can be of any

diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index ad6cd8eab6541cff50e46ddcf2d0e6fe125de961..de6c7f1e2a43cae4de6640d9c8e8cd7c1a944de9 100644 (file)
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -243,7 +243,7 @@ ptarray_to_GEOSCoordSeq(const POINTARRAY* pa, uint8_t fix_ring)
        uint32_t dims = 2;
        uint32_t i;
        int append_points = 0;
-       const POINT3DZ* p3d = NULL;
+       const POINT3D *p3d = NULL;
        const POINT2D* p2d = NULL;
        GEOSCoordSeq sq;
 
@@ -273,7 +273,7 @@ ptarray_to_GEOSCoordSeq(const POINTARRAY* pa, uint8_t fix_ring)
        {
                if (dims == 3)
                {
-                       p3d = getPoint3dz_cp(pa, i);
+                       p3d = getPoint3d_cp(pa, i);
                        p2d = (const POINT2D*)p3d;
                        LWDEBUGF(4, "Point: %g,%g,%g", p3d->x, p3d->y, p3d->z);
                }
@@ -293,7 +293,7 @@ ptarray_to_GEOSCoordSeq(const POINTARRAY* pa, uint8_t fix_ring)
        {
                if (dims == 3)
                {
-                       p3d = getPoint3dz_cp(pa, 0);
+                       p3d = getPoint3d_cp(pa, 0);
                        p2d = (const POINT2D*)p3d;
                }
                else

diff --git a/liblwgeom/lwinline.h b/liblwgeom/lwinline.h
index 35b4f0e5f9e633e50962b43750c065ef4d2beb36..880b51dc9d2c5abf3e2908b3efce2c7599b9b5c7 100644 (file)
--- a/liblwgeom/lwinline.h
+++ b/liblwgeom/lwinline.h
@@ -90,12 +90,33 @@ getPoint_internal(const POINTARRAY *pa, uint32_t n)
 static inline const POINT2D *
 getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
 {
-       if (!pa)
-               return 0;
-
        return (const POINT2D *)getPoint_internal(pa, n);
 }
 
+/**
+ * Returns a POINT2D pointer into the POINTARRAY serialized_ptlist,
+ * suitable for reading from. This is very high performance
+ * and declared const because you aren't allowed to muck with the
+ * values, only read them.
+ */
+static inline const POINT3D *
+getPoint3d_cp(const POINTARRAY *pa, uint32_t n)
+{
+       return (const POINT3D *)getPoint_internal(pa, n);
+}
+
+/**
+ * Returns a POINT2D pointer into the POINTARRAY serialized_ptlist,
+ * suitable for reading from. This is very high performance
+ * and declared const because you aren't allowed to muck with the
+ * values, only read them.
+ */
+static inline const POINT4D *
+getPoint4d_cp(const POINTARRAY *pa, uint32_t n)
+{
+       return (const POINT4D *)getPoint_internal(pa, n);
+}
+
 static inline LWPOINT *
 lwgeom_as_lwpoint(const LWGEOM *lwgeom)
 {

diff --git a/liblwgeom/lwout_geojson.c b/liblwgeom/lwout_geojson.c
index 36f14a295f38c7cd08cab69162405c7df3462ffb..f571279d3b1119eaf4255e9422a7fa438ff16f0b 100644 (file)
--- a/liblwgeom/lwout_geojson.c
+++ b/liblwgeom/lwout_geojson.c
@@ -752,8 +752,7 @@ pointArray_to_geojson(POINTARRAY *pa, char *output, int precision)
        {
                for (i=0; i<pa->npoints; i++)
                {
-                       const POINT3DZ *pt;
-                       pt = getPoint3dz_cp(pa, i);
+                       const POINT3D *pt = getPoint3d_cp(pa, i);
 
                        lwprint_double(
                            pt->x, precision, x, OUT_DOUBLE_BUFFER_SIZE);

diff --git a/liblwgeom/lwout_gml.c b/liblwgeom/lwout_gml.c
index f738f4dd887ed651c50fe2c8e136c402e0720092..d55d6dc19622173aa7fe8eca1f9da95e8a9ba54d 100644 (file)
--- a/liblwgeom/lwout_gml.c
+++ b/liblwgeom/lwout_gml.c
@@ -689,8 +689,7 @@ pointArray_toGML2(POINTARRAY *pa, char *output, int precision)
        {
                for (i=0; i<pa->npoints; i++)
                {
-                       const POINT3DZ *pt;
-                       pt = getPoint3dz_cp(pa, i);
+                       const POINT3D *pt = getPoint3d_cp(pa, i);
                        lwprint_double(
                            pt->x, precision, x, OUT_DOUBLE_BUFFER_SIZE);
                        lwprint_double(
@@ -1919,8 +1918,7 @@ pointArray_toGML3(POINTARRAY *pa, char *output, int precision, int opts)
        {
                for (i=0; i<pa->npoints; i++)
                {
-                       const POINT3DZ *pt;
-                       pt = getPoint3dz_cp(pa, i);
+                       const POINT3D *pt = getPoint3d_cp(pa, i);
 
                        lwprint_double(
                            pt->x, precision, x, OUT_DOUBLE_BUFFER_SIZE);