1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2018 Darafei Praliaskouski <me@komzpa.net>
22 * Copyright 2017-2018 Daniel Baston <dbaston@gmail.com>
23 * Copyright 2011 Sandro Santilli <strk@kbt.io>
24 * Copyright 2011 Paul Ramsey <pramsey@cleverelephant.ca>
25 * Copyright 2007-2008 Mark Cave-Ayland
26 * Copyright 2001-2006 Refractions Research Inc.
27 *
28 **********************************************************************/
29
30 #if PARANOIA_LEVEL > 0
31 #include <assert.h>
32 #endif
33
34 inline static double
distance2d_sqr_pt_pt(const POINT2D * p1,const POINT2D * p2)35 distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
36 {
37 double hside = p2->x - p1->x;
38 double vside = p2->y - p1->y;
39
40 return hside * hside + vside * vside;
41 }
42
43 inline static double
distance3d_sqr_pt_pt(const POINT3D * p1,const POINT3D * p2)44 distance3d_sqr_pt_pt(const POINT3D *p1, const POINT3D *p2)
45 {
46 double hside = p2->x - p1->x;
47 double vside = p2->y - p1->y;
48 double zside = p2->z - p1->z;
49
50 return hside * hside + vside * vside + zside * zside;
51 }
52
53 /*
54 * Size of point represeneted in the POINTARRAY
55 * 16 for 2d, 24 for 3d, 32 for 4d
56 */
57 static inline size_t
ptarray_point_size(const POINTARRAY * pa)58 ptarray_point_size(const POINTARRAY *pa)
59 {
60 return sizeof(double) * FLAGS_NDIMS(pa->flags);
61 }
62
63 /*
64 * Get a pointer to Nth point of a POINTARRAY
65 * You'll need to cast it to appropriate dimensioned point.
66 * Note that if you cast to a higher dimensional point you'll
67 * possibly corrupt the POINTARRAY.
68 *
69 * Casting to returned pointer to POINT2D* should be safe,
70 * as gserialized format always keeps the POINTARRAY pointer
71 * aligned to double boundary.
72 *
73 * WARNING: Don't cast this to a POINT!
74 * it would not be reliable due to memory alignment constraints
75 */
76 static inline uint8_t *
getPoint_internal(const POINTARRAY * pa,uint32_t n)77 getPoint_internal(const POINTARRAY *pa, uint32_t n)
78 {
79 size_t size;
80 uint8_t *ptr;
81
82 #if PARANOIA_LEVEL > 0
83 assert(pa);
84 assert(n <= pa->npoints);
85 assert(n <= pa->maxpoints);
86 #endif
87
88 size = ptarray_point_size(pa);
89 ptr = pa->serialized_pointlist + size * n;
90
91 return ptr;
92 }
93
94 /**
95 * Returns a POINT2D pointer into the POINTARRAY serialized_ptlist,
96 * suitable for reading from. This is very high performance
97 * and declared const because you aren't allowed to muck with the
98 * values, only read them.
99 */
100 static inline const POINT2D *
getPoint2d_cp(const POINTARRAY * pa,uint32_t n)101 getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
102 {
103 return (const POINT2D *)getPoint_internal(pa, n);
104 }
105
106 /**
107 * Returns a POINT3D pointer into the POINTARRAY serialized_ptlist,
108 * suitable for reading from. This is very high performance
109 * and declared const because you aren't allowed to muck with the
110 * values, only read them.
111 */
112 static inline const POINT3D *
getPoint3d_cp(const POINTARRAY * pa,uint32_t n)113 getPoint3d_cp(const POINTARRAY *pa, uint32_t n)
114 {
115 return (const POINT3D *)getPoint_internal(pa, n);
116 }
117
118 /**
119 * Returns a POINT4D pointer into the POINTARRAY serialized_ptlist,
120 * suitable for reading from. This is very high performance
121 * and declared const because you aren't allowed to muck with the
122 * values, only read them.
123 */
124 static inline const POINT4D *
getPoint4d_cp(const POINTARRAY * pa,uint32_t n)125 getPoint4d_cp(const POINTARRAY *pa, uint32_t n)
126 {
127 return (const POINT4D *)getPoint_internal(pa, n);
128 }
129
130 static inline LWPOINT *
lwgeom_as_lwpoint(const LWGEOM * lwgeom)131 lwgeom_as_lwpoint(const LWGEOM *lwgeom)
132 {
133 if (!lwgeom)
134 return NULL;
135 if (lwgeom->type == POINTTYPE)
136 return (LWPOINT *)lwgeom;
137 else
138 return NULL;
139 }
140
141 /**
142 * Return LWTYPE number
143 */
144 static inline uint32_t
lwgeom_get_type(const LWGEOM * geom)145 lwgeom_get_type(const LWGEOM *geom)
146 {
147 if (!geom)
148 return 0;
149 return geom->type;
150 }
151
152 static inline int
lwpoint_is_empty(const LWPOINT * point)153 lwpoint_is_empty(const LWPOINT *point)
154 {
155 return !point->point || point->point->npoints < 1;
156 }
157
158 static inline int
lwline_is_empty(const LWLINE * line)159 lwline_is_empty(const LWLINE *line)
160 {
161 return !line->points || line->points->npoints < 1;
162 }
163
164 static inline int
lwcircstring_is_empty(const LWCIRCSTRING * circ)165 lwcircstring_is_empty(const LWCIRCSTRING *circ)
166 {
167 return !circ->points || circ->points->npoints < 1;
168 }
169
170 static inline int
lwpoly_is_empty(const LWPOLY * poly)171 lwpoly_is_empty(const LWPOLY *poly)
172 {
173 return poly->nrings < 1 || !poly->rings || !poly->rings[0] || poly->rings[0]->npoints < 1;
174 }
175
176 static inline int
lwtriangle_is_empty(const LWTRIANGLE * triangle)177 lwtriangle_is_empty(const LWTRIANGLE *triangle)
178 {
179 return !triangle->points || triangle->points->npoints < 1;
180 }
181
182 static inline int lwgeom_is_empty(const LWGEOM *geom);
183
184 static inline int
lwcollection_is_empty(const LWCOLLECTION * col)185 lwcollection_is_empty(const LWCOLLECTION *col)
186 {
187 uint32_t i;
188 if (col->ngeoms == 0 || !col->geoms)
189 return LW_TRUE;
190 for (i = 0; i < col->ngeoms; i++)
191 {
192 if (!lwgeom_is_empty(col->geoms[i]))
193 return LW_FALSE;
194 }
195 return LW_TRUE;
196 }
197
198 /**
199 * Return true or false depending on whether a geometry is an "empty"
200 * geometry (no vertices members)
201 */
202 static inline int
lwgeom_is_empty(const LWGEOM * geom)203 lwgeom_is_empty(const LWGEOM *geom)
204 {
205 switch (geom->type)
206 {
207 case POINTTYPE:
208 return lwpoint_is_empty((LWPOINT *)geom);
209 break;
210 case LINETYPE:
211 return lwline_is_empty((LWLINE *)geom);
212 break;
213 case CIRCSTRINGTYPE:
214 return lwcircstring_is_empty((LWCIRCSTRING *)geom);
215 break;
216 case POLYGONTYPE:
217 return lwpoly_is_empty((LWPOLY *)geom);
218 break;
219 case TRIANGLETYPE:
220 return lwtriangle_is_empty((LWTRIANGLE *)geom);
221 break;
222 case MULTIPOINTTYPE:
223 case MULTILINETYPE:
224 case MULTIPOLYGONTYPE:
225 case COMPOUNDTYPE:
226 case CURVEPOLYTYPE:
227 case MULTICURVETYPE:
228 case MULTISURFACETYPE:
229 case POLYHEDRALSURFACETYPE:
230 case TINTYPE:
231 case COLLECTIONTYPE:
232 return lwcollection_is_empty((LWCOLLECTION *)geom);
233 break;
234 default:
235 return LW_FALSE;
236 break;
237 }
238 }
239
240 inline static uint64_t
uint64_interleave_2(uint64_t x,uint64_t y)241 uint64_interleave_2(uint64_t x, uint64_t y)
242 {
243 x = (x | (x << 16)) & 0x0000FFFF0000FFFFULL;
244 x = (x | (x << 8)) & 0x00FF00FF00FF00FFULL;
245 x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0FULL;
246 x = (x | (x << 2)) & 0x3333333333333333ULL;
247 x = (x | (x << 1)) & 0x5555555555555555ULL;
248
249 y = (y | (y << 16)) & 0x0000FFFF0000FFFFULL;
250 y = (y | (y << 8)) & 0x00FF00FF00FF00FFULL;
251 y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0FULL;
252 y = (y | (y << 2)) & 0x3333333333333333ULL;
253 y = (y | (y << 1)) & 0x5555555555555555ULL;
254
255 return x | (y << 1);
256 }
257
258 /* Based on https://github.com/rawrunprotected/hilbert_curves Public Domain code */
259 inline static uint64_t
uint32_hilbert(uint32_t px,uint32_t py)260 uint32_hilbert(uint32_t px, uint32_t py)
261 {
262 uint64_t x = px;
263 uint64_t y = py;
264
265 uint64_t A, B, C, D;
266 uint64_t a, b, c, d;
267 uint64_t i0, i1;
268
269 // Initial prefix scan round, prime with x and y
270 {
271 a = x ^ y;
272 b = 0xFFFFFFFFULL ^ a;
273 c = 0xFFFFFFFFULL ^ (x | y);
274 d = x & (y ^ 0xFFFFFFFFULL);
275
276 A = a | (b >> 1);
277 B = (a >> 1) ^ a;
278 C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
279 D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
280 }
281
282 {
283 a = A;
284 b = B;
285 c = C;
286 d = D;
287
288 A = ((a & (a >> 2)) ^ (b & (b >> 2)));
289 B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
290 C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
291 D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
292 }
293
294 {
295 a = A;
296 b = B;
297 c = C;
298 d = D;
299
300 A = ((a & (a >> 4)) ^ (b & (b >> 4)));
301 B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
302 C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
303 D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
304 }
305
306 {
307 a = A;
308 b = B;
309 c = C;
310 d = D;
311
312 A = ((a & (a >> 8)) ^ (b & (b >> 8)));
313 B = ((a & (b >> 8)) ^ (b & ((a ^ b) >> 8)));
314 C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
315 D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
316 }
317
318 {
319 a = A;
320 b = B;
321 c = C;
322 d = D;
323
324 C ^= ((a & (c >> 16)) ^ (b & (d >> 16)));
325 D ^= ((b & (c >> 16)) ^ ((a ^ b) & (d >> 16)));
326 }
327
328 // Undo transformation prefix scan
329 a = C ^ (C >> 1);
330 b = D ^ (D >> 1);
331
332 // Recover index bits
333 i0 = x ^ y;
334 i1 = b | (0xFFFFFFFFULL ^ (i0 | a));
335
336 return uint64_interleave_2(i0, i1);
337 }
338
339 /*
340 * This macro is based on PG_FREE_IF_COPY, except that it accepts two pointers.
341 * See PG_FREE_IF_COPY comment in src/include/fmgr.h in postgres source code
342 * for more details.
343 */
344 #define POSTGIS_FREE_IF_COPY_P(ptrsrc, ptrori) \
345 do \
346 { \
347 if ((Pointer)(ptrsrc) != (Pointer)(ptrori)) \
348 pfree(ptrsrc); \
349 } while (0)
350