/********************************************************************** * * 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 . * ********************************************************************** * * Copyright 2018 Darafei Praliaskouski * Copyright 2017-2018 Daniel Baston * Copyright 2011 Sandro Santilli * Copyright 2011 Paul Ramsey * Copyright 2007-2008 Mark Cave-Ayland * Copyright 2001-2006 Refractions Research Inc. * **********************************************************************/ #if PARANOIA_LEVEL > 0 #include #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; } /** * 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 POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n) { return (const POINT2D *)getPoint_internal(pa, n); } /** * Returns a POINT3D 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 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. */ 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; } /** * Return LWTYPE number */ 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; } /** * Return true or false depending on whether a geometry is an "empty" * geometry (no vertices members) */ 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)