lwinline.h

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 2018 Darafei Praliaskouski <me@komzpa.net>
 * Copyright 2017-2018 Daniel Baston <dbaston@gmail.com>
 * Copyright 2011 Sandro Santilli <strk@kbt.io>
 * Copyright 2011 Paul Ramsey <pramsey@cleverelephant.ca>
 * Copyright 2007-2008 Mark Cave-Ayland
 * Copyright 2001-2006 Refractions Research Inc.
 *
 **********************************************************************/
 
#if PARANOIA_LEVEL > 0
#include <assert.h>
#endif
 
inline static double
distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
{
        double hside = p2->x - p1->x;
        double vside = p2->y - p1->y;
 
        return hside * hside + vside * vside;
}
 
inline static double
distance3d_sqr_pt_pt(const POINT3D *p1, const POINT3D *p2)
{
        double hside = p2->x - p1->x;
        double vside = p2->y - p1->y;
        double zside = p2->z - p1->z;
 
        return hside * hside + vside * vside + zside * zside;
}
 
/*
 * Size of point represeneted in the POINTARRAY
 * 16 for 2d, 24 for 3d, 32 for 4d
 */
static inline size_t
ptarray_point_size(const POINTARRAY *pa)
{
        return sizeof(double) * FLAGS_NDIMS(pa->flags);
}
 
/*
 * Get a pointer to Nth point of a POINTARRAY
 * You'll need to cast it to appropriate dimensioned point.
 * Note that if you cast to a higher dimensional point you'll
 * possibly corrupt the POINTARRAY.
 *
 * Casting to returned pointer to POINT2D* should be safe,
 * as gserialized format always keeps the POINTARRAY pointer
 * aligned to double boundary.
 *
 * WARNING: Don't cast this to a POINT!
 * it would not be reliable due to memory alignment constraints
 */
static inline uint8_t *
getPoint_internal(const POINTARRAY *pa, uint32_t n)
{
        size_t size;
        uint8_t *ptr;
 
#if PARANOIA_LEVEL > 0
        assert(pa);
        assert(n <= pa->npoints);
        assert(n <= pa->maxpoints);
#endif
 
        size = ptarray_point_size(pa);
        ptr = pa->serialized_pointlist + size * n;
 
        return ptr;
}
 
static inline const POINT2D *
getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
{
        return (const POINT2D *)getPoint_internal(pa, n);
}
 
static inline const POINT3D *
getPoint3d_cp(const POINTARRAY *pa, uint32_t n)
{
        return (const POINT3D *)getPoint_internal(pa, n);
}
 
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)
{
        if (!lwgeom)
                return NULL;
        if (lwgeom->type == POINTTYPE)
                return (LWPOINT *)lwgeom;
        else
                return NULL;
}
 
static inline uint32_t
lwgeom_get_type(const LWGEOM *geom)
{
        if (!geom)
                return 0;
        return geom->type;
}
 
static inline int
lwpoint_is_empty(const LWPOINT *point)
{
        return !point->point || point->point->npoints < 1;
}
 
static inline int
lwline_is_empty(const LWLINE *line)
{
        return !line->points || line->points->npoints < 1;
}
 
static inline int
lwcircstring_is_empty(const LWCIRCSTRING *circ)
{
        return !circ->points || circ->points->npoints < 1;
}
 
static inline int
lwpoly_is_empty(const LWPOLY *poly)
{
        return poly->nrings < 1 || !poly->rings || !poly->rings[0] || poly->rings[0]->npoints < 1;
}
 
static inline int
lwtriangle_is_empty(const LWTRIANGLE *triangle)
{
        return !triangle->points || triangle->points->npoints < 1;
}
 
static inline int lwgeom_is_empty(const LWGEOM *geom);
 
static inline int
lwcollection_is_empty(const LWCOLLECTION *col)
{
        uint32_t i;
        if (col->ngeoms == 0 || !col->geoms)
                return LW_TRUE;
        for (i = 0; i < col->ngeoms; i++)
        {
                if (!lwgeom_is_empty(col->geoms[i]))
                        return LW_FALSE;
        }
        return LW_TRUE;
}
 
static inline int
lwgeom_is_empty(const LWGEOM *geom)
{
        switch (geom->type)
        {
        case POINTTYPE:
                return lwpoint_is_empty((LWPOINT *)geom);
                break;
        case LINETYPE:
                return lwline_is_empty((LWLINE *)geom);
                break;
        case CIRCSTRINGTYPE:
                return lwcircstring_is_empty((LWCIRCSTRING *)geom);
                break;
        case POLYGONTYPE:
                return lwpoly_is_empty((LWPOLY *)geom);
                break;
        case TRIANGLETYPE:
                return lwtriangle_is_empty((LWTRIANGLE *)geom);
                break;
        case MULTIPOINTTYPE:
        case MULTILINETYPE:
        case MULTIPOLYGONTYPE:
        case COMPOUNDTYPE:
        case CURVEPOLYTYPE:
        case MULTICURVETYPE:
        case MULTISURFACETYPE:
        case POLYHEDRALSURFACETYPE:
        case TINTYPE:
        case COLLECTIONTYPE:
                return lwcollection_is_empty((LWCOLLECTION *)geom);
                break;
        default:
                return LW_FALSE;
                break;
        }
}
 
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;
        uint64_t a, b, c, d;
        uint64_t i0, i1;
 
        // Initial prefix scan round, prime with x and y
        {
                a = x ^ y;
                b = 0xFFFFFFFFULL ^ a;
                c = 0xFFFFFFFFULL ^ (x | y);
                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;
        }
 
        {
                a = A;
                b = B;
                c = C;
                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)));
        }
 
        {
                a = A;
                b = B;
                c = C;
                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)));
        }
 
        {
                a = A;
                b = B;
                c = C;
                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)));
        }
 
        {
                a = A;
                b = B;
                c = C;
                d = D;
 
                C ^= ((a & (c >> 16)) ^ (b & (d >> 16)));
                D ^= ((b & (c >> 16)) ^ ((a ^ b) & (d >> 16)));
        }
 
        // Undo transformation prefix scan
        a = C ^ (C >> 1);
        b = D ^ (D >> 1);
 
        // Recover index bits
        i0 = x ^ y;
        i1 = b | (0xFFFFFFFFULL ^ (i0 | a));
 
        return uint64_interleave_2(i0, i1);
}
 
/*
 * This macro is based on PG_FREE_IF_COPY, except that it accepts two pointers.
 * See PG_FREE_IF_COPY comment in src/include/fmgr.h in postgres source code
 * for more details.
 */
#define POSTGIS_FREE_IF_COPY_P(ptrsrc, ptrori) \
        do \
        { \
                if ((Pointer)(ptrsrc) != (Pointer)(ptrori)) \
                        pfree(ptrsrc); \
        } while (0)