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 (C) 2012 Sandro Santilli <strk@kbt.io>
22 * Copyright (C) 2001-2006 Refractions Research Inc.
23 *
24 **********************************************************************/
25
26
27 /* basic LWPOLY manipulation */
28
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <math.h>
33 #include "liblwgeom_internal.h"
34 #include "lwgeom_log.h"
35
36
37 #define CHECK_POLY_RINGS_ZM 1
38
39 /* construct a new LWPOLY. arrays (points/points per ring) will NOT be copied
40 * use SRID=SRID_UNKNOWN for unknown SRID (will have 8bit type's S = 0)
41 */
42 LWPOLY*
lwpoly_construct(int srid,GBOX * bbox,uint32_t nrings,POINTARRAY ** points)43 lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
44 {
45 LWPOLY *result;
46 int hasz, hasm;
47 #ifdef CHECK_POLY_RINGS_ZM
48 char zm;
49 uint32_t i;
50 #endif
51
52 if ( nrings < 1 ) lwerror("lwpoly_construct: need at least 1 ring");
53
54 hasz = FLAGS_GET_Z(points[0]->flags);
55 hasm = FLAGS_GET_M(points[0]->flags);
56
57 #ifdef CHECK_POLY_RINGS_ZM
58 zm = FLAGS_GET_ZM(points[0]->flags);
59 for (i=1; i<nrings; i++)
60 {
61 if ( zm != FLAGS_GET_ZM(points[i]->flags) )
62 lwerror("lwpoly_construct: mixed dimensioned rings");
63 }
64 #endif
65
66 result = (LWPOLY*) lwalloc(sizeof(LWPOLY));
67 result->type = POLYGONTYPE;
68 result->flags = gflags(hasz, hasm, 0);
69 FLAGS_SET_BBOX(result->flags, bbox?1:0);
70 result->srid = srid;
71 result->nrings = nrings;
72 result->maxrings = nrings;
73 result->rings = points;
74 result->bbox = bbox;
75
76 return result;
77 }
78
79 LWPOLY*
lwpoly_construct_rectangle(char hasz,char hasm,POINT4D * p1,POINT4D * p2,POINT4D * p3,POINT4D * p4)80 lwpoly_construct_rectangle(char hasz, char hasm, POINT4D *p1, POINT4D *p2,
81 POINT4D *p3, POINT4D *p4)
82 {
83 POINTARRAY *pa = ptarray_construct_empty(hasz, hasm, 5);
84 LWPOLY *lwpoly = lwpoly_construct_empty(SRID_UNKNOWN, hasz, hasm);
85
86 ptarray_append_point(pa, p1, LW_TRUE);
87 ptarray_append_point(pa, p2, LW_TRUE);
88 ptarray_append_point(pa, p3, LW_TRUE);
89 ptarray_append_point(pa, p4, LW_TRUE);
90 ptarray_append_point(pa, p1, LW_TRUE);
91
92 lwpoly_add_ring(lwpoly, pa);
93
94 return lwpoly;
95 }
96
97 LWPOLY *
lwpoly_construct_envelope(int srid,double x1,double y1,double x2,double y2)98 lwpoly_construct_envelope(int srid, double x1, double y1, double x2, double y2)
99 {
100 POINT4D p1, p2, p3, p4;
101 LWPOLY *poly;
102
103 p1.x = x1;
104 p1.y = y1;
105 p2.x = x1;
106 p2.y = y2;
107 p3.x = x2;
108 p3.y = y2;
109 p4.x = x2;
110 p4.y = y1;
111
112 poly = lwpoly_construct_rectangle(0, 0, &p1, &p2, &p3, &p4);
113 lwgeom_set_srid(lwpoly_as_lwgeom(poly), srid);
114 lwgeom_add_bbox(lwpoly_as_lwgeom(poly));
115
116 return poly;
117 }
118
119 LWPOLY*
lwpoly_construct_circle(int srid,double x,double y,double radius,uint32_t segments_per_quarter,char exterior)120 lwpoly_construct_circle(int srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior)
121 {
122 const uint32_t segments = 4*segments_per_quarter;
123 double theta;
124 LWPOLY *lwpoly;
125 POINTARRAY *pa;
126 POINT4D pt;
127 uint32_t i;
128
129 if (segments_per_quarter == 0)
130 {
131 lwerror("Need at least one segment per quarter-circle.");
132 return NULL;
133 }
134
135 if (radius < 0)
136 {
137 lwerror("Radius must be positive.");
138 return NULL;
139 }
140
141 theta = 2*M_PI / segments;
142
143 lwpoly = lwpoly_construct_empty(srid, LW_FALSE, LW_FALSE);
144 pa = ptarray_construct_empty(LW_FALSE, LW_FALSE, segments + 1);
145
146 if (exterior)
147 radius *= sqrt(1 + pow(tan(theta/2), 2));
148
149 for (i = 0; i <= segments; i++)
150 {
151 pt.x = x + radius*sin(i * theta);
152 pt.y = y + radius*cos(i * theta);
153 ptarray_append_point(pa, &pt, LW_TRUE);
154 }
155
156 lwpoly_add_ring(lwpoly, pa);
157 return lwpoly;
158 }
159
160 LWPOLY*
lwpoly_construct_empty(int srid,char hasz,char hasm)161 lwpoly_construct_empty(int srid, char hasz, char hasm)
162 {
163 LWPOLY *result = lwalloc(sizeof(LWPOLY));
164 result->type = POLYGONTYPE;
165 result->flags = gflags(hasz,hasm,0);
166 result->srid = srid;
167 result->nrings = 0;
168 result->maxrings = 1; /* Allocate room for ring, just in case. */
169 result->rings = lwalloc(result->maxrings * sizeof(POINTARRAY*));
170 result->bbox = NULL;
171 return result;
172 }
173
174 void
lwpoly_free(LWPOLY * poly)175 lwpoly_free(LWPOLY* poly)
176 {
177 uint32_t t;
178
179 if (!poly) return;
180
181 if (poly->bbox) lwfree(poly->bbox);
182
183 if ( poly->rings )
184 {
185 for (t = 0; t < poly->nrings; t++)
186 if (poly->rings[t]) ptarray_free(poly->rings[t]);
187 lwfree(poly->rings);
188 }
189
190 lwfree(poly);
191 }
192
printLWPOLY(LWPOLY * poly)193 void printLWPOLY(LWPOLY *poly)
194 {
195 uint32_t t;
196 lwnotice("LWPOLY {");
197 lwnotice(" ndims = %i", (int)FLAGS_NDIMS(poly->flags));
198 lwnotice(" SRID = %i", (int)poly->srid);
199 lwnotice(" nrings = %i", (int)poly->nrings);
200 for (t=0; t<poly->nrings; t++)
201 {
202 lwnotice(" RING # %i :",t);
203 printPA(poly->rings[t]);
204 }
205 lwnotice("}");
206 }
207
208 /* @brief Clone LWLINE object. Serialized point lists are not copied.
209 *
210 * @see ptarray_clone
211 */
212 LWPOLY *
lwpoly_clone(const LWPOLY * g)213 lwpoly_clone(const LWPOLY *g)
214 {
215 uint32_t i;
216 LWPOLY *ret = lwalloc(sizeof(LWPOLY));
217 memcpy(ret, g, sizeof(LWPOLY));
218 ret->rings = lwalloc(sizeof(POINTARRAY *)*g->nrings);
219 for ( i = 0; i < g->nrings; i++ ) {
220 ret->rings[i] = ptarray_clone(g->rings[i]);
221 }
222 if ( g->bbox ) ret->bbox = gbox_copy(g->bbox);
223 return ret;
224 }
225
226 /* Deep clone LWPOLY object. POINTARRAY are copied, as is ring array */
227 LWPOLY *
lwpoly_clone_deep(const LWPOLY * g)228 lwpoly_clone_deep(const LWPOLY *g)
229 {
230 uint32_t i;
231 LWPOLY *ret = lwalloc(sizeof(LWPOLY));
232 memcpy(ret, g, sizeof(LWPOLY));
233 if ( g->bbox ) ret->bbox = gbox_copy(g->bbox);
234 ret->rings = lwalloc(sizeof(POINTARRAY *)*g->nrings);
235 for ( i = 0; i < ret->nrings; i++ )
236 {
237 ret->rings[i] = ptarray_clone_deep(g->rings[i]);
238 }
239 FLAGS_SET_READONLY(ret->flags,0);
240 return ret;
241 }
242
243 /**
244 * Add a ring to a polygon. Point array will be referenced, not copied.
245 */
246 int
lwpoly_add_ring(LWPOLY * poly,POINTARRAY * pa)247 lwpoly_add_ring(LWPOLY *poly, POINTARRAY *pa)
248 {
249 if( ! poly || ! pa )
250 return LW_FAILURE;
251
252 /* We have used up our storage, add some more. */
253 if( poly->nrings >= poly->maxrings )
254 {
255 int new_maxrings = 2 * (poly->nrings + 1);
256 poly->rings = lwrealloc(poly->rings, new_maxrings * sizeof(POINTARRAY*));
257 poly->maxrings = new_maxrings;
258 }
259
260 /* Add the new ring entry. */
261 poly->rings[poly->nrings] = pa;
262 poly->nrings++;
263
264 return LW_SUCCESS;
265 }
266
267 void
lwpoly_force_clockwise(LWPOLY * poly)268 lwpoly_force_clockwise(LWPOLY *poly)
269 {
270 uint32_t i;
271
272 /* No-op empties */
273 if ( lwpoly_is_empty(poly) )
274 return;
275
276 /* External ring */
277 if ( ptarray_isccw(poly->rings[0]) )
278 ptarray_reverse_in_place(poly->rings[0]);
279
280 /* Internal rings */
281 for (i=1; i<poly->nrings; i++)
282 if ( ! ptarray_isccw(poly->rings[i]) )
283 ptarray_reverse_in_place(poly->rings[i]);
284
285 }
286
287 int
lwpoly_is_clockwise(LWPOLY * poly)288 lwpoly_is_clockwise(LWPOLY *poly)
289 {
290 uint32_t i;
291
292 if ( lwpoly_is_empty(poly) )
293 return LW_TRUE;
294
295 if ( ptarray_isccw(poly->rings[0]) )
296 return LW_FALSE;
297
298 for ( i = 1; i < poly->nrings; i++)
299 if ( !ptarray_isccw(poly->rings[i]) )
300 return LW_FALSE;
301
302 return LW_TRUE;
303 }
304
305 void
lwpoly_release(LWPOLY * lwpoly)306 lwpoly_release(LWPOLY *lwpoly)
307 {
308 lwgeom_release(lwpoly_as_lwgeom(lwpoly));
309 }
310
311 LWPOLY *
lwpoly_segmentize2d(const LWPOLY * poly,double dist)312 lwpoly_segmentize2d(const LWPOLY *poly, double dist)
313 {
314 POINTARRAY **newrings;
315 uint32_t i;
316
317 newrings = lwalloc(sizeof(POINTARRAY *)*poly->nrings);
318 for (i=0; i<poly->nrings; i++)
319 {
320 newrings[i] = ptarray_segmentize2d(poly->rings[i], dist);
321 if ( ! newrings[i] )
322 {
323 uint32_t j = 0;
324 for (j = 0; j < i; j++)
325 ptarray_free(newrings[j]);
326 lwfree(newrings);
327 return NULL;
328 }
329 }
330 return lwpoly_construct(poly->srid, NULL,
331 poly->nrings, newrings);
332 }
333
334 /*
335 * check coordinate equality
336 * ring and coordinate order is considered
337 */
338 char
lwpoly_same(const LWPOLY * p1,const LWPOLY * p2)339 lwpoly_same(const LWPOLY *p1, const LWPOLY *p2)
340 {
341 uint32_t i;
342
343 if ( p1->nrings != p2->nrings ) return 0;
344 for (i=0; i<p1->nrings; i++)
345 {
346 if ( ! ptarray_same(p1->rings[i], p2->rings[i]) )
347 return 0;
348 }
349 return 1;
350 }
351
352 /*
353 * Construct a polygon from a LWLINE being
354 * the shell and an array of LWLINE (possibly NULL) being holes.
355 * Pointarrays from intput geoms are cloned.
356 * SRID must be the same for each input line.
357 * Input lines must have at least 4 points, and be closed.
358 */
359 LWPOLY *
lwpoly_from_lwlines(const LWLINE * shell,uint32_t nholes,const LWLINE ** holes)360 lwpoly_from_lwlines(const LWLINE *shell,
361 uint32_t nholes, const LWLINE **holes)
362 {
363 uint32_t nrings;
364 POINTARRAY **rings = lwalloc((nholes+1)*sizeof(POINTARRAY *));
365 int srid = shell->srid;
366 LWPOLY *ret;
367
368 if ( shell->points->npoints < 4 )
369 lwerror("lwpoly_from_lwlines: shell must have at least 4 points");
370 if ( ! ptarray_is_closed_2d(shell->points) )
371 lwerror("lwpoly_from_lwlines: shell must be closed");
372 rings[0] = ptarray_clone_deep(shell->points);
373
374 for (nrings=1; nrings<=nholes; nrings++)
375 {
376 const LWLINE *hole = holes[nrings-1];
377
378 if ( hole->srid != srid )
379 lwerror("lwpoly_from_lwlines: mixed SRIDs in input lines");
380
381 if ( hole->points->npoints < 4 )
382 lwerror("lwpoly_from_lwlines: holes must have at least 4 points");
383 if ( ! ptarray_is_closed_2d(hole->points) )
384 lwerror("lwpoly_from_lwlines: holes must be closed");
385
386 rings[nrings] = ptarray_clone_deep(hole->points);
387 }
388
389 ret = lwpoly_construct(srid, NULL, nrings, rings);
390 return ret;
391 }
392
393 LWPOLY*
lwpoly_force_dims(const LWPOLY * poly,int hasz,int hasm)394 lwpoly_force_dims(const LWPOLY *poly, int hasz, int hasm)
395 {
396 LWPOLY *polyout;
397
398 /* Return 2D empty */
399 if( lwpoly_is_empty(poly) )
400 {
401 polyout = lwpoly_construct_empty(poly->srid, hasz, hasm);
402 }
403 else
404 {
405 POINTARRAY **rings = NULL;
406 uint32_t i;
407 rings = lwalloc(sizeof(POINTARRAY*) * poly->nrings);
408 for( i = 0; i < poly->nrings; i++ )
409 {
410 rings[i] = ptarray_force_dims(poly->rings[i], hasz, hasm);
411 }
412 polyout = lwpoly_construct(poly->srid, NULL, poly->nrings, rings);
413 }
414 polyout->type = poly->type;
415 return polyout;
416 }
417
lwpoly_is_empty(const LWPOLY * poly)418 int lwpoly_is_empty(const LWPOLY *poly)
419 {
420 if ( (poly->nrings < 1) || (!poly->rings) || (!poly->rings[0]) || (poly->rings[0]->npoints < 1) )
421 return LW_TRUE;
422 return LW_FALSE;
423 }
424
lwpoly_count_vertices(LWPOLY * poly)425 uint32_t lwpoly_count_vertices(LWPOLY *poly)
426 {
427 uint32_t i = 0;
428 uint32_t v = 0; /* vertices */
429 assert(poly);
430 for ( i = 0; i < poly->nrings; i ++ )
431 {
432 v += poly->rings[i]->npoints;
433 }
434 return v;
435 }
436
437 /**
438 * Find the area of the outer ring - sum (area of inner rings).
439 */
440 double
lwpoly_area(const LWPOLY * poly)441 lwpoly_area(const LWPOLY *poly)
442 {
443 double poly_area = 0.0;
444 uint32_t i;
445
446 if ( ! poly )
447 lwerror("lwpoly_area called with null polygon pointer!");
448
449 for ( i=0; i < poly->nrings; i++ )
450 {
451 POINTARRAY *ring = poly->rings[i];
452 double ringarea = 0.0;
453
454 /* Empty or messed-up ring. */
455 if ( ring->npoints < 3 )
456 continue;
457
458 ringarea = fabs(ptarray_signed_area(ring));
459 if ( i == 0 ) /* Outer ring, positive area! */
460 poly_area += ringarea;
461 else /* Inner ring, negative area! */
462 poly_area -= ringarea;
463 }
464
465 return poly_area;
466 }
467
468
469 /**
470 * Compute the sum of polygon rings length.
471 * Could use a more numerically stable calculator...
472 */
473 double
lwpoly_perimeter(const LWPOLY * poly)474 lwpoly_perimeter(const LWPOLY *poly)
475 {
476 double result=0.0;
477 uint32_t i;
478
479 LWDEBUGF(2, "in lwgeom_polygon_perimeter (%d rings)", poly->nrings);
480
481 for (i=0; i<poly->nrings; i++)
482 result += ptarray_length(poly->rings[i]);
483
484 return result;
485 }
486
487 /**
488 * Compute the sum of polygon rings length (forcing 2d computation).
489 * Could use a more numerically stable calculator...
490 */
491 double
lwpoly_perimeter_2d(const LWPOLY * poly)492 lwpoly_perimeter_2d(const LWPOLY *poly)
493 {
494 double result=0.0;
495 uint32_t i;
496
497 LWDEBUGF(2, "in lwgeom_polygon_perimeter (%d rings)", poly->nrings);
498
499 for (i=0; i<poly->nrings; i++)
500 result += ptarray_length_2d(poly->rings[i]);
501
502 return result;
503 }
504
505 int
lwpoly_is_closed(const LWPOLY * poly)506 lwpoly_is_closed(const LWPOLY *poly)
507 {
508 uint32_t i = 0;
509
510 if ( poly->nrings == 0 )
511 return LW_TRUE;
512
513 for ( i = 0; i < poly->nrings; i++ )
514 {
515 if (FLAGS_GET_Z(poly->flags))
516 {
517 if ( ! ptarray_is_closed_3d(poly->rings[i]) )
518 return LW_FALSE;
519 }
520 else
521 {
522 if ( ! ptarray_is_closed_2d(poly->rings[i]) )
523 return LW_FALSE;
524 }
525 }
526
527 return LW_TRUE;
528 }
529
530 int
lwpoly_startpoint(const LWPOLY * poly,POINT4D * pt)531 lwpoly_startpoint(const LWPOLY* poly, POINT4D* pt)
532 {
533 if ( poly->nrings < 1 )
534 return LW_FAILURE;
535 return ptarray_startpoint(poly->rings[0], pt);
536 }
537
538 int
lwpoly_contains_point(const LWPOLY * poly,const POINT2D * pt)539 lwpoly_contains_point(const LWPOLY *poly, const POINT2D *pt)
540 {
541 uint32_t i;
542
543 if ( lwpoly_is_empty(poly) )
544 return LW_FALSE;
545
546 if ( ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE )
547 return LW_FALSE;
548
549 for ( i = 1; i < poly->nrings; i++ )
550 {
551 if ( ptarray_contains_point(poly->rings[i], pt) == LW_INSIDE )
552 return LW_FALSE;
553 }
554 return LW_TRUE;
555 }
556
557
558