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) 2001-2006 Refractions Research Inc.
22 * Copyright (C) 2017 Sandro Santilli <strk@kbt.io>
23 * Copyright (C) 2018 Daniel Baston <dbaston@gmail.com>
24 *
25 **********************************************************************/
26
27
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <stdarg.h>
31 #include <string.h>
32
33 #include "../postgis_config.h"
34
35 /*#define POSTGIS_DEBUG_LEVEL 3*/
36
37 #include "lwgeom_log.h"
38
39 #include "liblwgeom_internal.h"
40
41 LWGEOM *pta_unstroke(const POINTARRAY *points, int32_t srid);
42 LWGEOM* lwline_unstroke(const LWLINE *line);
43 LWGEOM* lwpolygon_unstroke(const LWPOLY *poly);
44 LWGEOM* lwmline_unstroke(const LWMLINE *mline);
45 LWGEOM* lwmpolygon_unstroke(const LWMPOLY *mpoly);
46 LWGEOM* lwcollection_unstroke(const LWCOLLECTION *c);
47 LWGEOM* lwgeom_unstroke(const LWGEOM *geom);
48
49
50 /*
51 * Determines (recursively in the case of collections) whether the geometry
52 * contains at least on arc geometry or segment.
53 */
54 int
lwgeom_has_arc(const LWGEOM * geom)55 lwgeom_has_arc(const LWGEOM *geom)
56 {
57 LWCOLLECTION *col;
58 uint32_t i;
59
60 LWDEBUG(2, "lwgeom_has_arc called.");
61
62 switch (geom->type)
63 {
64 case POINTTYPE:
65 case LINETYPE:
66 case POLYGONTYPE:
67 case TRIANGLETYPE:
68 case MULTIPOINTTYPE:
69 case MULTILINETYPE:
70 case MULTIPOLYGONTYPE:
71 case POLYHEDRALSURFACETYPE:
72 case TINTYPE:
73 return LW_FALSE;
74 case CIRCSTRINGTYPE:
75 return LW_TRUE;
76 /* It's a collection that MAY contain an arc */
77 default:
78 col = (LWCOLLECTION *)geom;
79 for (i=0; i<col->ngeoms; i++)
80 {
81 if (lwgeom_has_arc(col->geoms[i]) == LW_TRUE)
82 return LW_TRUE;
83 }
84 return LW_FALSE;
85 }
86 }
87
88 int
lwgeom_type_arc(const LWGEOM * geom)89 lwgeom_type_arc(const LWGEOM *geom)
90 {
91 switch (geom->type)
92 {
93 case COMPOUNDTYPE:
94 case CIRCSTRINGTYPE:
95 case CURVEPOLYTYPE:
96 case MULTISURFACETYPE:
97 case MULTICURVETYPE:
98 return LW_TRUE;
99 default:
100 return LW_FALSE;
101 }
102 }
103
104 /*******************************************************************************
105 * Begin curve segmentize functions
106 ******************************************************************************/
107
interpolate_arc(double angle,double a1,double a2,double a3,double zm1,double zm2,double zm3)108 static double interpolate_arc(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
109 {
110 LWDEBUGF(4,"angle %.05g a1 %.05g a2 %.05g a3 %.05g zm1 %.05g zm2 %.05g zm3 %.05g",angle,a1,a2,a3,zm1,zm2,zm3);
111 /* Counter-clockwise sweep */
112 if ( a1 < a2 )
113 {
114 if ( angle <= a2 )
115 return zm1 + (zm2-zm1) * (angle-a1) / (a2-a1);
116 else
117 return zm2 + (zm3-zm2) * (angle-a2) / (a3-a2);
118 }
119 /* Clockwise sweep */
120 else
121 {
122 if ( angle >= a2 )
123 return zm1 + (zm2-zm1) * (a1-angle) / (a1-a2);
124 else
125 return zm2 + (zm3-zm2) * (a2-angle) / (a2-a3);
126 }
127 }
128
129 /* Compute the angle covered by a single segment such that
130 * a given number of segments per quadrant is achieved. */
angle_increment_using_segments_per_quad(double tol)131 static double angle_increment_using_segments_per_quad(double tol)
132 {
133 double increment;
134 int perQuad = rint(tol);
135 // error out if tol != perQuad ? (not-round)
136 if ( perQuad != tol )
137 {
138 lwerror("lwarc_linearize: segments per quadrant must be an integer value, got %.15g", tol, perQuad);
139 return -1;
140 }
141 if ( perQuad < 1 )
142 {
143 lwerror("lwarc_linearize: segments per quadrant must be at least 1, got %d", perQuad);
144 return -1;
145 }
146 increment = fabs(M_PI_2 / perQuad);
147 LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g (%g degrees)", perQuad, increment, increment*180/M_PI);
148
149 return increment;
150 }
151
152 /* Compute the angle covered by a single quadrant such that
153 * the segment deviates from the arc by no more than a given
154 * amount. */
angle_increment_using_max_deviation(double max_deviation,double radius)155 static double angle_increment_using_max_deviation(double max_deviation, double radius)
156 {
157 double increment, halfAngle, maxErr;
158 if ( max_deviation <= 0 )
159 {
160 lwerror("lwarc_linearize: max deviation must be bigger than 0, got %.15g", max_deviation);
161 return -1;
162 }
163
164 /*
165 * Ref: https://en.wikipedia.org/wiki/Sagitta_(geometry)
166 *
167 * An arc "sagitta" (distance between middle point of arc and
168 * middle point of corresponding chord) is defined as:
169 *
170 * sagitta = radius * ( 1 - cos( angle ) );
171 *
172 * We want our sagitta to be at most "tolerance" long,
173 * and we want to find out angle, so we use the inverse
174 * formula:
175 *
176 * tol = radius * ( 1 - cos( angle ) );
177 * 1 - cos( angle ) = tol/radius
178 * - cos( angle ) = tol/radius - 1
179 * cos( angle ) = - tol/radius + 1
180 * angle = acos( 1 - tol/radius )
181 *
182 * Constraints: 1.0 - tol/radius must be between -1 and 1
183 * which means tol must be between 0 and 2 times
184 * the radius, which makes sense as you cannot have a
185 * sagitta bigger than twice the radius!
186 *
187 */
188 maxErr = max_deviation;
189 if ( maxErr > radius * 2 )
190 {
191 maxErr = radius * 2;
192 LWDEBUGF(2,
193 "lwarc_linearize: tolerance %g is too big, "
194 "using arc-max 2 * radius == %g",
195 max_deviation,
196 maxErr);
197 }
198 do {
199 halfAngle = acos( 1.0 - maxErr / radius );
200 /* TODO: avoid a loop here, going rather straight to
201 * a minimum angle value */
202 if ( halfAngle != 0 ) break;
203 LWDEBUGF(2, "lwarc_linearize: tolerance %g is too small for this arc"
204 " to compute approximation angle, doubling it", maxErr);
205 maxErr *= 2;
206 } while(1);
207 increment = 2 * halfAngle;
208 LWDEBUGF(2,
209 "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)",
210 max_deviation,
211 radius,
212 halfAngle,
213 increment,
214 increment * 180 / M_PI);
215
216 return increment;
217 }
218
219 /* Check that a given angle is positive and, if so, take
220 * it to be the angle covered by a single segment. */
angle_increment_using_max_angle(double tol)221 static double angle_increment_using_max_angle(double tol)
222 {
223 if ( tol <= 0 )
224 {
225 lwerror("lwarc_linearize: max angle must be bigger than 0, got %.15g", tol);
226 return -1;
227 }
228
229 return tol;
230 }
231
232
233 /**
234 * Segmentize an arc
235 *
236 * Does not add the final vertex
237 *
238 * @param to POINTARRAY to append segmentized vertices to
239 * @param p1 first point defining the arc
240 * @param p2 second point defining the arc
241 * @param p3 third point defining the arc
242 * @param tol tolerance, semantic driven by tolerance_type
243 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
244 * @param flags LW_LINEARIZE_FLAGS
245 *
246 * @return number of points appended (0 if collinear),
247 * or -1 on error (lwerror would be called).
248 *
249 */
250 static int
lwarc_linearize(POINTARRAY * to,const POINT4D * p1,const POINT4D * p2,const POINT4D * p3,double tol,LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,int flags)251 lwarc_linearize(POINTARRAY *to,
252 const POINT4D *p1, const POINT4D *p2, const POINT4D *p3,
253 double tol, LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
254 int flags)
255 {
256 POINT2D center;
257 POINT2D *t1 = (POINT2D*)p1;
258 POINT2D *t2 = (POINT2D*)p2;
259 POINT2D *t3 = (POINT2D*)p3;
260 POINT4D pt;
261 int p2_side = 0;
262 int clockwise = LW_TRUE;
263 double radius; /* Arc radius */
264 double increment; /* Angle per segment */
265 double angle_shift = 0;
266 double a1, a2, a3;
267 POINTARRAY *pa;
268 int is_circle = LW_FALSE;
269 int points_added = 0;
270 int reverse = 0;
271 int segments = 0;
272
273 LWDEBUG(2, "lwarc_linearize called.");
274
275 LWDEBUGF(2, " curve is CIRCULARSTRING(%.15g %.15f, %.15f %.15f, %.15f %15f)",
276 t1->x, t1->y, t2->x, t2->y, t3->x, t3->y);
277
278 p2_side = lw_segment_side(t1, t3, t2);
279
280 LWDEBUGF(2, " p2 side is %d", p2_side);
281
282 /* Force counterclockwise scan if SYMMETRIC operation is requested */
283 if ( p2_side == -1 && flags & LW_LINEARIZE_FLAG_SYMMETRIC )
284 {
285 /* swap p1-p3 */
286 t1 = (POINT2D*)p3;
287 t3 = (POINT2D*)p1;
288 p1 = (POINT4D*)t1;
289 p3 = (POINT4D*)t3;
290 p2_side = 1;
291 reverse = 1;
292 }
293
294 radius = lw_arc_center(t1, t2, t3, ¢er);
295 LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
296
297 /* Matched start/end points imply circle */
298 if ( p1->x == p3->x && p1->y == p3->y )
299 is_circle = LW_TRUE;
300
301 /* Negative radius signals straight line, p1/p2/p3 are collinear */
302 if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
303 return 0;
304
305 /* The side of the p1/p3 line that p2 falls on dictates the sweep
306 direction from p1 to p3. */
307 if ( p2_side == -1 )
308 clockwise = LW_TRUE;
309 else
310 clockwise = LW_FALSE;
311
312 /* Compute the increment (angle per segment) depending on
313 * our tolerance type. */
314 switch(tolerance_type)
315 {
316 case LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD:
317 increment = angle_increment_using_segments_per_quad(tol);
318 break;
319 case LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION:
320 increment = angle_increment_using_max_deviation(tol, radius);
321 break;
322 case LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE:
323 increment = angle_increment_using_max_angle(tol);
324 break;
325 default:
326 lwerror("lwarc_linearize: unsupported tolerance type %d", tolerance_type);
327 return -1;
328 }
329
330 if (increment < 0)
331 {
332 /* Error occurred in increment calculation somewhere
333 * (lwerror already called)
334 */
335 return -1;
336 }
337
338 /* Angles of each point that defines the arc section */
339 a1 = atan2(p1->y - center.y, p1->x - center.x);
340 a2 = atan2(p2->y - center.y, p2->x - center.x);
341 a3 = atan2(p3->y - center.y, p3->x - center.x);
342
343 LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
344 a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);
345
346 /* Calculate total arc angle, in radians */
347 double total_angle = clockwise ? a1 - a3 : a3 - a1;
348 if ( total_angle <= 0 ) total_angle += M_PI * 2;
349
350 /* At extreme tolerance values (very low or very high, depending on
351 * the semantic) we may cause our arc to collapse. In this case,
352 * we want shrink the increment enough so that we get two segments
353 * for a standard arc, or three segments for a complete circle. */
354 int min_segs = is_circle ? 3 : 2;
355 segments = ceil(total_angle / increment);
356 if (segments < min_segs)
357 {
358 segments = min_segs;
359 increment = total_angle / min_segs;
360 }
361
362 if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
363 {
364 LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg", total_angle * 180 / M_PI);
365
366 if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
367 {
368 /* Number of complete steps */
369 segments = trunc(total_angle / increment);
370
371 /* Figure out the angle remainder, i.e. the amount of the angle
372 * that is left after we can take no more complete angle
373 * increments. */
374 double angle_remainder = total_angle - (increment * segments);
375
376 /* Shift the starting angle by half of the remainder. This
377 * will have the effect of evenly distributing the remainder
378 * among the first and last segments in the arc. */
379 angle_shift = angle_remainder / 2.0;
380
381 LWDEBUGF(2,
382 "lwarc_linearize RETAIN_ANGLE operation requested - "
383 "total angle %g, steps %d, increment %g, remainder %g",
384 total_angle * 180 / M_PI,
385 segments,
386 increment * 180 / M_PI,
387 angle_remainder * 180 / M_PI);
388 }
389 else
390 {
391 /* Number of segments in output */
392 segments = ceil(total_angle / increment);
393 /* Tweak increment to be regular for all the arc */
394 increment = total_angle / segments;
395
396 LWDEBUGF(2,
397 "lwarc_linearize SYMMETRIC operation requested - "
398 "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d - I:%g",
399 total_angle * 180 / M_PI,
400 p1->x,
401 p1->y,
402 center.x,
403 center.y,
404 p3->x,
405 p3->y,
406 segments,
407 increment * 180 / M_PI);
408 }
409 }
410
411 /* p2 on left side => clockwise sweep */
412 if ( clockwise )
413 {
414 LWDEBUG(2, " Clockwise sweep");
415 increment *= -1;
416 angle_shift *= -1;
417 /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
418 if ( a3 > a1 )
419 a3 -= 2.0 * M_PI;
420 if ( a2 > a1 )
421 a2 -= 2.0 * M_PI;
422 }
423 /* p2 on right side => counter-clockwise sweep */
424 else
425 {
426 LWDEBUG(2, " Counterclockwise sweep");
427 /* Adjust a3 up so we can increment from a1 to a3 cleanly */
428 if ( a3 < a1 )
429 a3 += 2.0 * M_PI;
430 if ( a2 < a1 )
431 a2 += 2.0 * M_PI;
432 }
433
434 /* Override angles for circle case */
435 if( is_circle )
436 {
437 increment = fabs(increment);
438 segments = ceil(total_angle / increment);
439 if (segments < 3)
440 {
441 segments = 3;
442 increment = total_angle / 3;
443 }
444 a3 = a1 + 2.0 * M_PI;
445 a2 = a1 + M_PI;
446 clockwise = LW_FALSE;
447 angle_shift = 0.0;
448 }
449
450 LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
451 angle_shift * 180/M_PI, increment * 180/M_PI);
452
453 if ( reverse )
454 {
455 /* Append points in order to a temporary POINTARRAY and
456 * reverse them before writing to the output POINTARRAY. */
457 const int capacity = 8; /* TODO: compute exactly ? */
458 pa = ptarray_construct_empty(ptarray_has_z(to), ptarray_has_m(to), capacity);
459 }
460 else
461 {
462 /* Append points directly to the output POINTARRAY,
463 * starting with p1. */
464 pa = to;
465
466 ptarray_append_point(pa, p1, LW_FALSE);
467 ++points_added;
468 }
469
470 /* Sweep from a1 to a3 */
471 int seg_start = 1; /* First point is added manually */
472 int seg_end = segments;
473 if (angle_shift != 0.0)
474 {
475 /* When we have extra angles we need to add the extra segments at the
476 * start and end that cover those parts of the arc */
477 seg_start = 0;
478 seg_end = segments + 1;
479 }
480 LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
481 a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
482 for (int s = seg_start; s < seg_end; s++)
483 {
484 double angle = a1 + increment * s + angle_shift;
485 LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
486 pt.x = center.x + radius * cos(angle);
487 pt.y = center.y + radius * sin(angle);
488 pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
489 pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
490 ptarray_append_point(pa, &pt, LW_FALSE);
491 ++points_added;
492 }
493
494 /* Ensure the final point is EXACTLY the same as the first for the circular case */
495 if ( is_circle )
496 {
497 ptarray_remove_point(pa, pa->npoints - 1);
498 ptarray_append_point(pa, p1, LW_FALSE);
499 }
500
501 if ( reverse )
502 {
503 int i;
504 ptarray_append_point(to, p3, LW_FALSE);
505 for ( i=pa->npoints; i>0; i-- ) {
506 getPoint4d_p(pa, i-1, &pt);
507 ptarray_append_point(to, &pt, LW_FALSE);
508 }
509 ptarray_free(pa);
510 }
511
512 return points_added;
513 }
514
515 /*
516 * @param icurve input curve
517 * @param tol tolerance, semantic driven by tolerance_type
518 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
519 * @param flags see flags in lwarc_linearize
520 *
521 * @return a newly allocated LWLINE
522 */
523 static LWLINE *
lwcircstring_linearize(const LWCIRCSTRING * icurve,double tol,LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,int flags)524 lwcircstring_linearize(const LWCIRCSTRING *icurve, double tol,
525 LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
526 int flags)
527 {
528 LWLINE *oline;
529 POINTARRAY *ptarray;
530 uint32_t i, j;
531 POINT4D p1, p2, p3, p4;
532 int ret;
533
534 LWDEBUGF(2, "lwcircstring_linearize called., dim = %d", icurve->points->flags);
535
536 ptarray = ptarray_construct_empty(FLAGS_GET_Z(icurve->points->flags), FLAGS_GET_M(icurve->points->flags), 64);
537
538 for (i = 2; i < icurve->points->npoints; i+=2)
539 {
540 LWDEBUGF(3, "lwcircstring_linearize: arc ending at point %d", i);
541
542 getPoint4d_p(icurve->points, i - 2, &p1);
543 getPoint4d_p(icurve->points, i - 1, &p2);
544 getPoint4d_p(icurve->points, i, &p3);
545
546 ret = lwarc_linearize(ptarray, &p1, &p2, &p3, tol, tolerance_type, flags);
547 if ( ret > 0 )
548 {
549 LWDEBUGF(3, "lwcircstring_linearize: generated %d points", ptarray->npoints);
550 }
551 else if ( ret == 0 )
552 {
553 LWDEBUG(3, "lwcircstring_linearize: points are colinear, returning curve points as line");
554
555 for (j = i - 2 ; j < i ; j++)
556 {
557 getPoint4d_p(icurve->points, j, &p4);
558 ptarray_append_point(ptarray, &p4, LW_TRUE);
559 }
560 }
561 else
562 {
563 /* An error occurred, lwerror should have been called by now */
564 ptarray_free(ptarray);
565 return NULL;
566 }
567 }
568 getPoint4d_p(icurve->points, icurve->points->npoints-1, &p1);
569 ptarray_append_point(ptarray, &p1, LW_FALSE);
570
571 oline = lwline_construct(icurve->srid, NULL, ptarray);
572 return oline;
573 }
574
575 /*
576 * @param icompound input compound curve
577 * @param tol tolerance, semantic driven by tolerance_type
578 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
579 * @param flags see flags in lwarc_linearize
580 *
581 * @return a newly allocated LWLINE
582 */
583 static LWLINE *
lwcompound_linearize(const LWCOMPOUND * icompound,double tol,LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,int flags)584 lwcompound_linearize(const LWCOMPOUND *icompound, double tol,
585 LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
586 int flags)
587 {
588 LWGEOM *geom;
589 POINTARRAY *ptarray = NULL;
590 LWLINE *tmp = NULL;
591 uint32_t i, j;
592 POINT4D p;
593
594 LWDEBUG(2, "lwcompound_stroke called.");
595
596 ptarray = ptarray_construct_empty(FLAGS_GET_Z(icompound->flags), FLAGS_GET_M(icompound->flags), 64);
597
598 for (i = 0; i < icompound->ngeoms; i++)
599 {
600 geom = icompound->geoms[i];
601 if (geom->type == CIRCSTRINGTYPE)
602 {
603 tmp = lwcircstring_linearize((LWCIRCSTRING *)geom, tol, tolerance_type, flags);
604 for (j = 0; j < tmp->points->npoints; j++)
605 {
606 getPoint4d_p(tmp->points, j, &p);
607 ptarray_append_point(ptarray, &p, LW_TRUE);
608 }
609 lwline_free(tmp);
610 }
611 else if (geom->type == LINETYPE)
612 {
613 tmp = (LWLINE *)geom;
614 for (j = 0; j < tmp->points->npoints; j++)
615 {
616 getPoint4d_p(tmp->points, j, &p);
617 ptarray_append_point(ptarray, &p, LW_TRUE);
618 }
619 }
620 else
621 {
622 lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(geom->type));
623 return NULL;
624 }
625 }
626
627 ptarray_remove_repeated_points_in_place(ptarray, 0.0, 2);
628 return lwline_construct(icompound->srid, NULL, ptarray);
629 }
630
631
632 /*
633 * @param icompound input curve polygon
634 * @param tol tolerance, semantic driven by tolerance_type
635 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
636 * @param flags see flags in lwarc_linearize
637 *
638 * @return a newly allocated LWPOLY
639 */
640 static LWPOLY *
lwcurvepoly_linearize(const LWCURVEPOLY * curvepoly,double tol,LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,int flags)641 lwcurvepoly_linearize(const LWCURVEPOLY *curvepoly, double tol,
642 LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
643 int flags)
644 {
645 LWPOLY *ogeom;
646 LWGEOM *tmp;
647 LWLINE *line;
648 POINTARRAY **ptarray;
649 uint32_t i;
650
651 LWDEBUG(2, "lwcurvepoly_linearize called.");
652
653 ptarray = lwalloc(sizeof(POINTARRAY *)*curvepoly->nrings);
654
655 for (i = 0; i < curvepoly->nrings; i++)
656 {
657 tmp = curvepoly->rings[i];
658 if (tmp->type == CIRCSTRINGTYPE)
659 {
660 line = lwcircstring_linearize((LWCIRCSTRING *)tmp, tol, tolerance_type, flags);
661 ptarray[i] = ptarray_clone_deep(line->points);
662 lwline_free(line);
663 }
664 else if (tmp->type == LINETYPE)
665 {
666 line = (LWLINE *)tmp;
667 ptarray[i] = ptarray_clone_deep(line->points);
668 }
669 else if (tmp->type == COMPOUNDTYPE)
670 {
671 line = lwcompound_linearize((LWCOMPOUND *)tmp, tol, tolerance_type, flags);
672 ptarray[i] = ptarray_clone_deep(line->points);
673 lwline_free(line);
674 }
675 else
676 {
677 lwerror("Invalid ring type found in CurvePoly.");
678 return NULL;
679 }
680 }
681
682 ogeom = lwpoly_construct(curvepoly->srid, NULL, curvepoly->nrings, ptarray);
683 return ogeom;
684 }
685
686 /* Kept for backward compatibility - TODO: drop */
687 LWPOLY *
lwcurvepoly_stroke(const LWCURVEPOLY * curvepoly,uint32_t perQuad)688 lwcurvepoly_stroke(const LWCURVEPOLY *curvepoly, uint32_t perQuad)
689 {
690 return lwcurvepoly_linearize(curvepoly, perQuad, LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD, 0);
691 }
692
693
694 /**
695 * @param mcurve input compound curve
696 * @param tol tolerance, semantic driven by tolerance_type
697 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
698 * @param flags see flags in lwarc_linearize
699 *
700 * @return a newly allocated LWMLINE
701 */
702 static LWMLINE *
lwmcurve_linearize(const LWMCURVE * mcurve,double tol,LW_LINEARIZE_TOLERANCE_TYPE type,int flags)703 lwmcurve_linearize(const LWMCURVE *mcurve, double tol,
704 LW_LINEARIZE_TOLERANCE_TYPE type,
705 int flags)
706 {
707 LWMLINE *ogeom;
708 LWGEOM **lines;
709 uint32_t i;
710
711 LWDEBUGF(2, "lwmcurve_linearize called, geoms=%d, dim=%d.", mcurve->ngeoms, FLAGS_NDIMS(mcurve->flags));
712
713 lines = lwalloc(sizeof(LWGEOM *)*mcurve->ngeoms);
714
715 for (i = 0; i < mcurve->ngeoms; i++)
716 {
717 const LWGEOM *tmp = mcurve->geoms[i];
718 if (tmp->type == CIRCSTRINGTYPE)
719 {
720 lines[i] = (LWGEOM *)lwcircstring_linearize((LWCIRCSTRING *)tmp, tol, type, flags);
721 }
722 else if (tmp->type == LINETYPE)
723 {
724 lines[i] = (LWGEOM *)lwline_construct(mcurve->srid, NULL, ptarray_clone_deep(((LWLINE *)tmp)->points));
725 }
726 else if (tmp->type == COMPOUNDTYPE)
727 {
728 lines[i] = (LWGEOM *)lwcompound_linearize((LWCOMPOUND *)tmp, tol, type, flags);
729 }
730 else
731 {
732 lwerror("Unsupported geometry found in MultiCurve.");
733 return NULL;
734 }
735 }
736
737 ogeom = (LWMLINE *)lwcollection_construct(MULTILINETYPE, mcurve->srid, NULL, mcurve->ngeoms, lines);
738 return ogeom;
739 }
740
741 /**
742 * @param msurface input multi surface
743 * @param tol tolerance, semantic driven by tolerance_type
744 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
745 * @param flags see flags in lwarc_linearize
746 *
747 * @return a newly allocated LWMPOLY
748 */
749 static LWMPOLY *
lwmsurface_linearize(const LWMSURFACE * msurface,double tol,LW_LINEARIZE_TOLERANCE_TYPE type,int flags)750 lwmsurface_linearize(const LWMSURFACE *msurface, double tol,
751 LW_LINEARIZE_TOLERANCE_TYPE type,
752 int flags)
753 {
754 LWMPOLY *ogeom;
755 LWGEOM *tmp;
756 LWPOLY *poly;
757 LWGEOM **polys;
758 POINTARRAY **ptarray;
759 uint32_t i, j;
760
761 LWDEBUG(2, "lwmsurface_linearize called.");
762
763 polys = lwalloc(sizeof(LWGEOM *)*msurface->ngeoms);
764
765 for (i = 0; i < msurface->ngeoms; i++)
766 {
767 tmp = msurface->geoms[i];
768 if (tmp->type == CURVEPOLYTYPE)
769 {
770 polys[i] = (LWGEOM *)lwcurvepoly_linearize((LWCURVEPOLY *)tmp, tol, type, flags);
771 }
772 else if (tmp->type == POLYGONTYPE)
773 {
774 poly = (LWPOLY *)tmp;
775 ptarray = lwalloc(sizeof(POINTARRAY *)*poly->nrings);
776 for (j = 0; j < poly->nrings; j++)
777 {
778 ptarray[j] = ptarray_clone_deep(poly->rings[j]);
779 }
780 polys[i] = (LWGEOM *)lwpoly_construct(msurface->srid, NULL, poly->nrings, ptarray);
781 }
782 }
783 ogeom = (LWMPOLY *)lwcollection_construct(MULTIPOLYGONTYPE, msurface->srid, NULL, msurface->ngeoms, polys);
784 return ogeom;
785 }
786
787 /**
788 * @param collection input geometry collection
789 * @param tol tolerance, semantic driven by tolerance_type
790 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
791 * @param flags see flags in lwarc_linearize
792 *
793 * @return a newly allocated LWCOLLECTION
794 */
795 static LWCOLLECTION *
lwcollection_linearize(const LWCOLLECTION * collection,double tol,LW_LINEARIZE_TOLERANCE_TYPE type,int flags)796 lwcollection_linearize(const LWCOLLECTION *collection, double tol,
797 LW_LINEARIZE_TOLERANCE_TYPE type,
798 int flags)
799 {
800 LWCOLLECTION *ocol;
801 LWGEOM *tmp;
802 LWGEOM **geoms;
803 uint32_t i;
804
805 LWDEBUG(2, "lwcollection_linearize called.");
806
807 geoms = lwalloc(sizeof(LWGEOM *)*collection->ngeoms);
808
809 for (i=0; i<collection->ngeoms; i++)
810 {
811 tmp = collection->geoms[i];
812 switch (tmp->type)
813 {
814 case CIRCSTRINGTYPE:
815 geoms[i] = (LWGEOM *)lwcircstring_linearize((LWCIRCSTRING *)tmp, tol, type, flags);
816 break;
817 case COMPOUNDTYPE:
818 geoms[i] = (LWGEOM *)lwcompound_linearize((LWCOMPOUND *)tmp, tol, type, flags);
819 break;
820 case CURVEPOLYTYPE:
821 geoms[i] = (LWGEOM *)lwcurvepoly_linearize((LWCURVEPOLY *)tmp, tol, type, flags);
822 break;
823 case MULTICURVETYPE:
824 case MULTISURFACETYPE:
825 case COLLECTIONTYPE:
826 geoms[i] = (LWGEOM *)lwcollection_linearize((LWCOLLECTION *)tmp, tol, type, flags);
827 break;
828 default:
829 geoms[i] = lwgeom_clone_deep(tmp);
830 break;
831 }
832 }
833 ocol = lwcollection_construct(COLLECTIONTYPE, collection->srid, NULL, collection->ngeoms, geoms);
834 return ocol;
835 }
836
837 LWGEOM *
lwcurve_linearize(const LWGEOM * geom,double tol,LW_LINEARIZE_TOLERANCE_TYPE type,int flags)838 lwcurve_linearize(const LWGEOM *geom, double tol,
839 LW_LINEARIZE_TOLERANCE_TYPE type,
840 int flags)
841 {
842 LWGEOM * ogeom = NULL;
843 switch (geom->type)
844 {
845 case CIRCSTRINGTYPE:
846 ogeom = (LWGEOM *)lwcircstring_linearize((LWCIRCSTRING *)geom, tol, type, flags);
847 break;
848 case COMPOUNDTYPE:
849 ogeom = (LWGEOM *)lwcompound_linearize((LWCOMPOUND *)geom, tol, type, flags);
850 break;
851 case CURVEPOLYTYPE:
852 ogeom = (LWGEOM *)lwcurvepoly_linearize((LWCURVEPOLY *)geom, tol, type, flags);
853 break;
854 case MULTICURVETYPE:
855 ogeom = (LWGEOM *)lwmcurve_linearize((LWMCURVE *)geom, tol, type, flags);
856 break;
857 case MULTISURFACETYPE:
858 ogeom = (LWGEOM *)lwmsurface_linearize((LWMSURFACE *)geom, tol, type, flags);
859 break;
860 case COLLECTIONTYPE:
861 ogeom = (LWGEOM *)lwcollection_linearize((LWCOLLECTION *)geom, tol, type, flags);
862 break;
863 default:
864 ogeom = lwgeom_clone_deep(geom);
865 }
866 return ogeom;
867 }
868
869 /* Kept for backward compatibility - TODO: drop */
870 LWGEOM *
lwgeom_stroke(const LWGEOM * geom,uint32_t perQuad)871 lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)
872 {
873 return lwcurve_linearize(geom, perQuad, LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD, 0);
874 }
875
876 /**
877 * Return ABC angle in radians
878 * TODO: move to lwalgorithm
879 */
880 static double
lw_arc_angle(const POINT2D * a,const POINT2D * b,const POINT2D * c)881 lw_arc_angle(const POINT2D *a, const POINT2D *b, const POINT2D *c)
882 {
883 POINT2D ab, cb;
884
885 ab.x = b->x - a->x;
886 ab.y = b->y - a->y;
887
888 cb.x = b->x - c->x;
889 cb.y = b->y - c->y;
890
891 double dot = (ab.x * cb.x + ab.y * cb.y); /* dot product */
892 double cross = (ab.x * cb.y - ab.y * cb.x); /* cross product */
893
894 double alpha = atan2(cross, dot);
895
896 return alpha;
897 }
898
899 /**
900 * Returns LW_TRUE if b is on the arc formed by a1/a2/a3, but not within
901 * that portion already described by a1/a2/a3
902 */
pt_continues_arc(const POINT4D * a1,const POINT4D * a2,const POINT4D * a3,const POINT4D * b)903 static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b)
904 {
905 POINT2D center;
906 POINT2D *t1 = (POINT2D*)a1;
907 POINT2D *t2 = (POINT2D*)a2;
908 POINT2D *t3 = (POINT2D*)a3;
909 POINT2D *tb = (POINT2D*)b;
910 double radius = lw_arc_center(t1, t2, t3, ¢er);
911 double b_distance, diff;
912
913 /* Co-linear a1/a2/a3 */
914 if ( radius < 0.0 )
915 return LW_FALSE;
916
917 b_distance = distance2d_pt_pt(tb, ¢er);
918 diff = fabs(radius - b_distance);
919 LWDEBUGF(4, "circle_radius=%g, b_distance=%g, diff=%g, percentage=%g", radius, b_distance, diff, diff/radius);
920
921 /* Is the point b on the circle? */
922 if ( diff < EPSILON_SQLMM )
923 {
924 int a2_side = lw_segment_side(t1, t3, t2);
925 int b_side = lw_segment_side(t1, t3, tb);
926 double angle1 = lw_arc_angle(t1, t2, t3);
927 double angle2 = lw_arc_angle(t2, t3, tb);
928
929 /* Is the angle similar to the previous one ? */
930 diff = fabs(angle1 - angle2);
931 LWDEBUGF(4, " angle1: %g, angle2: %g, diff:%g", angle1, angle2, diff);
932 if ( diff > EPSILON_SQLMM )
933 {
934 return LW_FALSE;
935 }
936
937 /* Is the point b on the same side of a1/a3 as the mid-point a2 is? */
938 /* If not, it's in the unbounded part of the circle, so it continues the arc, return true. */
939 if ( b_side != a2_side )
940 return LW_TRUE;
941 }
942 return LW_FALSE;
943 }
944
945 static LWGEOM *
linestring_from_pa(const POINTARRAY * pa,int32_t srid,int start,int end)946 linestring_from_pa(const POINTARRAY *pa, int32_t srid, int start, int end)
947 {
948 int i = 0, j = 0;
949 POINT4D p;
950 POINTARRAY *pao = ptarray_construct(ptarray_has_z(pa), ptarray_has_m(pa), end-start+2);
951 LWDEBUGF(4, "srid=%d, start=%d, end=%d", srid, start, end);
952 for( i = start; i < end + 2; i++ )
953 {
954 getPoint4d_p(pa, i, &p);
955 ptarray_set_point4d(pao, j++, &p);
956 }
957 return lwline_as_lwgeom(lwline_construct(srid, NULL, pao));
958 }
959
960 static LWGEOM *
circstring_from_pa(const POINTARRAY * pa,int32_t srid,int start,int end)961 circstring_from_pa(const POINTARRAY *pa, int32_t srid, int start, int end)
962 {
963
964 POINT4D p0, p1, p2;
965 POINTARRAY *pao = ptarray_construct(ptarray_has_z(pa), ptarray_has_m(pa), 3);
966 LWDEBUGF(4, "srid=%d, start=%d, end=%d", srid, start, end);
967 getPoint4d_p(pa, start, &p0);
968 ptarray_set_point4d(pao, 0, &p0);
969 getPoint4d_p(pa, (start+end+1)/2, &p1);
970 ptarray_set_point4d(pao, 1, &p1);
971 getPoint4d_p(pa, end+1, &p2);
972 ptarray_set_point4d(pao, 2, &p2);
973 return lwcircstring_as_lwgeom(lwcircstring_construct(srid, NULL, pao));
974 }
975
976 static LWGEOM *
geom_from_pa(const POINTARRAY * pa,int32_t srid,int is_arc,int start,int end)977 geom_from_pa(const POINTARRAY *pa, int32_t srid, int is_arc, int start, int end)
978 {
979 LWDEBUGF(4, "srid=%d, is_arc=%d, start=%d, end=%d", srid, is_arc, start, end);
980 if ( is_arc )
981 return circstring_from_pa(pa, srid, start, end);
982 else
983 return linestring_from_pa(pa, srid, start, end);
984 }
985
986 LWGEOM *
pta_unstroke(const POINTARRAY * points,int32_t srid)987 pta_unstroke(const POINTARRAY *points, int32_t srid)
988 {
989 int i = 0, j, k;
990 POINT4D a1, a2, a3, b;
991 POINT4D first, center;
992 char *edges_in_arcs;
993 int found_arc = LW_FALSE;
994 int current_arc = 1;
995 int num_edges;
996 int edge_type; /* non-zero if edge is part of an arc */
997 int start, end;
998 LWCOLLECTION *outcol;
999 /* Minimum number of edges, per quadrant, required to define an arc */
1000 const unsigned int min_quad_edges = 2;
1001
1002 /* Die on null input */
1003 if ( ! points )
1004 lwerror("pta_unstroke called with null pointarray");
1005
1006 /* Null on empty input? */
1007 if ( points->npoints == 0 )
1008 return NULL;
1009
1010 /* We can't desegmentize anything shorter than four points */
1011 if ( points->npoints < 4 )
1012 {
1013 /* Return a linestring here*/
1014 lwerror("pta_unstroke needs implementation for npoints < 4");
1015 }
1016
1017 /* Allocate our result array of vertices that are part of arcs */
1018 num_edges = points->npoints - 1;
1019 edges_in_arcs = lwalloc(num_edges + 1);
1020 memset(edges_in_arcs, 0, num_edges + 1);
1021
1022 /* We make a candidate arc of the first two edges, */
1023 /* And then see if the next edge follows it */
1024 while( i < num_edges-2 )
1025 {
1026 unsigned int arc_edges;
1027 double num_quadrants;
1028 double angle;
1029
1030 found_arc = LW_FALSE;
1031 /* Make candidate arc */
1032 getPoint4d_p(points, i , &a1);
1033 getPoint4d_p(points, i+1, &a2);
1034 getPoint4d_p(points, i+2, &a3);
1035 memcpy(&first, &a1, sizeof(POINT4D));
1036
1037 for( j = i+3; j < num_edges+1; j++ )
1038 {
1039 LWDEBUGF(4, "i=%d, j=%d", i, j);
1040 getPoint4d_p(points, j, &b);
1041 /* Does this point fall on our candidate arc? */
1042 if ( pt_continues_arc(&a1, &a2, &a3, &b) )
1043 {
1044 /* Yes. Mark this edge and the two preceding it as arc components */
1045 LWDEBUGF(4, "pt_continues_arc #%d", current_arc);
1046 found_arc = LW_TRUE;
1047 for ( k = j-1; k > j-4; k-- )
1048 edges_in_arcs[k] = current_arc;
1049 }
1050 else
1051 {
1052 /* No. So we're done with this candidate arc */
1053 LWDEBUG(4, "pt_continues_arc = false");
1054 current_arc++;
1055 break;
1056 }
1057
1058 memcpy(&a1, &a2, sizeof(POINT4D));
1059 memcpy(&a2, &a3, sizeof(POINT4D));
1060 memcpy(&a3, &b, sizeof(POINT4D));
1061 }
1062 /* Jump past all the edges that were added to the arc */
1063 if ( found_arc )
1064 {
1065 /* Check if an arc was composed by enough edges to be
1066 * really considered an arc
1067 * See http://trac.osgeo.org/postgis/ticket/2420
1068 */
1069 arc_edges = j - 1 - i;
1070 LWDEBUGF(4, "arc defined by %d edges found", arc_edges);
1071 if ( first.x == b.x && first.y == b.y ) {
1072 LWDEBUG(4, "arc is a circle");
1073 num_quadrants = 4;
1074 }
1075 else {
1076 lw_arc_center((POINT2D*)&first, (POINT2D*)&b, (POINT2D*)&a1, (POINT2D*)¢er);
1077 angle = lw_arc_angle((POINT2D*)&first, (POINT2D*)¢er, (POINT2D*)&b);
1078 int p2_side = lw_segment_side((POINT2D*)&first, (POINT2D*)&a1, (POINT2D*)&b);
1079 if ( p2_side >= 0 ) angle = -angle;
1080
1081 if ( angle < 0 ) angle = 2 * M_PI + angle;
1082 num_quadrants = ( 4 * angle ) / ( 2 * M_PI );
1083 LWDEBUGF(4, "arc angle (%g %g, %g %g, %g %g) is %g (side is %d), quadrants:%g", first.x, first.y, center.x, center.y, b.x, b.y, angle, p2_side, num_quadrants);
1084 }
1085 /* a1 is first point, b is last point */
1086 if ( arc_edges < min_quad_edges * num_quadrants ) {
1087 LWDEBUGF(4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants);
1088 for ( k = j-1; k >= i; k-- )
1089 edges_in_arcs[k] = 0;
1090 }
1091
1092 i = j-1;
1093 }
1094 else
1095 {
1096 /* Mark this edge as a linear edge */
1097 edges_in_arcs[i] = 0;
1098 i = i+1;
1099 }
1100 }
1101
1102 #if POSTGIS_DEBUG_LEVEL > 3
1103 {
1104 char *edgestr = lwalloc(num_edges+1);
1105 for ( i = 0; i < num_edges; i++ )
1106 {
1107 if ( edges_in_arcs[i] )
1108 edgestr[i] = 48 + edges_in_arcs[i];
1109 else
1110 edgestr[i] = '.';
1111 }
1112 edgestr[num_edges] = 0;
1113 LWDEBUGF(3, "edge pattern %s", edgestr);
1114 lwfree(edgestr);
1115 }
1116 #endif
1117
1118 start = 0;
1119 edge_type = edges_in_arcs[0];
1120 outcol = lwcollection_construct_empty(COMPOUNDTYPE, srid, ptarray_has_z(points), ptarray_has_m(points));
1121 for( i = 1; i < num_edges; i++ )
1122 {
1123 if( edge_type != edges_in_arcs[i] )
1124 {
1125 end = i - 1;
1126 lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1127 start = i;
1128 edge_type = edges_in_arcs[i];
1129 }
1130 }
1131 lwfree(edges_in_arcs); /* not needed anymore */
1132
1133 /* Roll out last item */
1134 end = num_edges - 1;
1135 lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
1136
1137 /* Strip down to singleton if only one entry */
1138 if ( outcol->ngeoms == 1 )
1139 {
1140 LWGEOM *outgeom = outcol->geoms[0];
1141 outcol->ngeoms = 0; lwcollection_free(outcol);
1142 return outgeom;
1143 }
1144 return lwcollection_as_lwgeom(outcol);
1145 }
1146
1147
1148 LWGEOM *
lwline_unstroke(const LWLINE * line)1149 lwline_unstroke(const LWLINE *line)
1150 {
1151 LWDEBUG(2, "lwline_unstroke called.");
1152
1153 if ( line->points->npoints < 4 ) return lwline_as_lwgeom(lwline_clone_deep(line));
1154 else return pta_unstroke(line->points, line->srid);
1155 }
1156
1157 LWGEOM *
lwpolygon_unstroke(const LWPOLY * poly)1158 lwpolygon_unstroke(const LWPOLY *poly)
1159 {
1160 LWGEOM **geoms;
1161 uint32_t i, hascurve = 0;
1162
1163 LWDEBUG(2, "lwpolygon_unstroke called.");
1164
1165 geoms = lwalloc(sizeof(LWGEOM *)*poly->nrings);
1166 for (i=0; i<poly->nrings; i++)
1167 {
1168 geoms[i] = pta_unstroke(poly->rings[i], poly->srid);
1169 if (geoms[i]->type == CIRCSTRINGTYPE || geoms[i]->type == COMPOUNDTYPE)
1170 {
1171 hascurve = 1;
1172 }
1173 }
1174 if (hascurve == 0)
1175 {
1176 for (i=0; i<poly->nrings; i++)
1177 {
1178 lwfree(geoms[i]); /* TODO: should this be lwgeom_free instead ? */
1179 }
1180 return lwgeom_clone_deep((LWGEOM *)poly);
1181 }
1182
1183 return (LWGEOM *)lwcollection_construct(CURVEPOLYTYPE, poly->srid, NULL, poly->nrings, geoms);
1184 }
1185
1186 LWGEOM *
lwmline_unstroke(const LWMLINE * mline)1187 lwmline_unstroke(const LWMLINE *mline)
1188 {
1189 LWGEOM **geoms;
1190 uint32_t i, hascurve = 0;
1191
1192 LWDEBUG(2, "lwmline_unstroke called.");
1193
1194 geoms = lwalloc(sizeof(LWGEOM *)*mline->ngeoms);
1195 for (i=0; i<mline->ngeoms; i++)
1196 {
1197 geoms[i] = lwline_unstroke((LWLINE *)mline->geoms[i]);
1198 if (geoms[i]->type == CIRCSTRINGTYPE || geoms[i]->type == COMPOUNDTYPE)
1199 {
1200 hascurve = 1;
1201 }
1202 }
1203 if (hascurve == 0)
1204 {
1205 for (i=0; i<mline->ngeoms; i++)
1206 {
1207 lwfree(geoms[i]); /* TODO: should this be lwgeom_free instead ? */
1208 }
1209 return lwgeom_clone_deep((LWGEOM *)mline);
1210 }
1211 return (LWGEOM *)lwcollection_construct(MULTICURVETYPE, mline->srid, NULL, mline->ngeoms, geoms);
1212 }
1213
1214 LWGEOM *
lwmpolygon_unstroke(const LWMPOLY * mpoly)1215 lwmpolygon_unstroke(const LWMPOLY *mpoly)
1216 {
1217 LWGEOM **geoms;
1218 uint32_t i, hascurve = 0;
1219
1220 LWDEBUG(2, "lwmpoly_unstroke called.");
1221
1222 geoms = lwalloc(sizeof(LWGEOM *)*mpoly->ngeoms);
1223 for (i=0; i<mpoly->ngeoms; i++)
1224 {
1225 geoms[i] = lwpolygon_unstroke((LWPOLY *)mpoly->geoms[i]);
1226 if (geoms[i]->type == CURVEPOLYTYPE)
1227 {
1228 hascurve = 1;
1229 }
1230 }
1231 if (hascurve == 0)
1232 {
1233 for (i=0; i<mpoly->ngeoms; i++)
1234 {
1235 lwfree(geoms[i]); /* TODO: should this be lwgeom_free instead ? */
1236 }
1237 return lwgeom_clone_deep((LWGEOM *)mpoly);
1238 }
1239 return (LWGEOM *)lwcollection_construct(MULTISURFACETYPE, mpoly->srid, NULL, mpoly->ngeoms, geoms);
1240 }
1241
1242 LWGEOM *
lwcollection_unstroke(const LWCOLLECTION * c)1243 lwcollection_unstroke(const LWCOLLECTION *c)
1244 {
1245 LWCOLLECTION *ret = lwalloc(sizeof(LWCOLLECTION));
1246 memcpy(ret, c, sizeof(LWCOLLECTION));
1247
1248 if (c->ngeoms > 0)
1249 {
1250 uint32_t i;
1251 ret->geoms = lwalloc(sizeof(LWGEOM *)*c->ngeoms);
1252 for (i=0; i < c->ngeoms; i++)
1253 {
1254 ret->geoms[i] = lwgeom_unstroke(c->geoms[i]);
1255 }
1256 if (c->bbox)
1257 {
1258 ret->bbox = gbox_copy(c->bbox);
1259 }
1260 }
1261 else
1262 {
1263 ret->bbox = NULL;
1264 ret->geoms = NULL;
1265 }
1266 return (LWGEOM *)ret;
1267 }
1268
1269
1270 LWGEOM *
lwgeom_unstroke(const LWGEOM * geom)1271 lwgeom_unstroke(const LWGEOM *geom)
1272 {
1273 LWDEBUG(2, "lwgeom_unstroke called.");
1274
1275 switch (geom->type)
1276 {
1277 case LINETYPE:
1278 return lwline_unstroke((LWLINE *)geom);
1279 case POLYGONTYPE:
1280 return lwpolygon_unstroke((LWPOLY *)geom);
1281 case MULTILINETYPE:
1282 return lwmline_unstroke((LWMLINE *)geom);
1283 case MULTIPOLYGONTYPE:
1284 return lwmpolygon_unstroke((LWMPOLY *)geom);
1285 case COLLECTIONTYPE:
1286 return lwcollection_unstroke((LWCOLLECTION *)geom);
1287 default:
1288 return lwgeom_clone_deep(geom);
1289 }
1290 }
1291
1292