PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwsegmentize.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: lwsegmentize.c 13134 2014-12-01 08:47:21Z strk $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  * Copyright 2001-2006 Refractions Research Inc.
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <stdarg.h>
16 #include <string.h>
17 
18 #include "liblwgeom_internal.h"
19 
20 /* #define POSTGIS_DEBUG_LEVEL 4 */
21 
22 #include "lwgeom_log.h"
23 
24 
25 LWMLINE *lwmcurve_segmentize(LWMCURVE *mcurve, uint32_t perQuad);
26 LWMPOLY *lwmsurface_segmentize(LWMSURFACE *msurface, uint32_t perQuad);
27 LWCOLLECTION *lwcollection_segmentize(LWCOLLECTION *collection, uint32_t perQuad);
28 
29 LWGEOM *pta_desegmentize(POINTARRAY *points, int type, int srid);
35 
36 
37 /*
38  * Determines (recursively in the case of collections) whether the geometry
39  * contains at least on arc geometry or segment.
40  */
41 int
42 lwgeom_has_arc(const LWGEOM *geom)
43 {
44  LWCOLLECTION *col;
45  int i;
46 
47  LWDEBUG(2, "lwgeom_has_arc called.");
48 
49  switch (geom->type)
50  {
51  case POINTTYPE:
52  case LINETYPE:
53  case POLYGONTYPE:
54  case TRIANGLETYPE:
55  case MULTIPOINTTYPE:
56  case MULTILINETYPE:
57  case MULTIPOLYGONTYPE:
59  case TINTYPE:
60  return LW_FALSE;
61  case CIRCSTRINGTYPE:
62  return LW_TRUE;
63  /* It's a collection that MAY contain an arc */
64  default:
65  col = (LWCOLLECTION *)geom;
66  for (i=0; i<col->ngeoms; i++)
67  {
68  if (lwgeom_has_arc(col->geoms[i]) == LW_TRUE)
69  return LW_TRUE;
70  }
71  return LW_FALSE;
72  }
73 }
74 
75 
76 
77 /*******************************************************************************
78  * Begin curve segmentize functions
79  ******************************************************************************/
80 
81 static double interpolate_arc(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
82 {
83  LWDEBUGF(4,"angle %.05g a1 %.05g a2 %.05g a3 %.05g zm1 %.05g zm2 %.05g zm3 %.05g",angle,a1,a2,a3,zm1,zm2,zm3);
84  /* Counter-clockwise sweep */
85  if ( a1 < a2 )
86  {
87  if ( angle <= a2 )
88  return zm1 + (zm2-zm1) * (angle-a1) / (a2-a1);
89  else
90  return zm2 + (zm3-zm2) * (angle-a2) / (a3-a2);
91  }
92  /* Clockwise sweep */
93  else
94  {
95  if ( angle >= a2 )
96  return zm1 + (zm2-zm1) * (a1-angle) / (a1-a2);
97  else
98  return zm2 + (zm3-zm2) * (a2-angle) / (a2-a3);
99  }
100 }
101 
102 static POINTARRAY *
103 lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad)
104 {
105  POINT2D center;
106  POINT2D *t1 = (POINT2D*)p1;
107  POINT2D *t2 = (POINT2D*)p2;
108  POINT2D *t3 = (POINT2D*)p3;
109  POINT4D pt;
110  int p2_side = 0;
111  int clockwise = LW_TRUE;
112  double radius; /* Arc radius */
113  double increment; /* Angle per segment */
114  double a1, a2, a3, angle;
115  POINTARRAY *pa;
116  int is_circle = LW_FALSE;
117 
118  LWDEBUG(2, "lwcircle_calculate_gbox called.");
119 
120  radius = lw_arc_center(t1, t2, t3, &center);
121  p2_side = lw_segment_side(t1, t3, t2);
122 
123  /* Matched start/end points imply circle */
124  if ( p1->x == p3->x && p1->y == p3->y )
125  is_circle = LW_TRUE;
126 
127  /* Negative radius signals straight line, p1/p2/p3 are colinear */
128  if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
129  return NULL;
130 
131  /* The side of the p1/p3 line that p2 falls on dictates the sweep
132  direction from p1 to p3. */
133  if ( p2_side == -1 )
134  clockwise = LW_TRUE;
135  else
136  clockwise = LW_FALSE;
137 
138  increment = fabs(M_PI_2 / perQuad);
139 
140  /* Angles of each point that defines the arc section */
141  a1 = atan2(p1->y - center.y, p1->x - center.x);
142  a2 = atan2(p2->y - center.y, p2->x - center.x);
143  a3 = atan2(p3->y - center.y, p3->x - center.x);
144 
145  /* p2 on left side => clockwise sweep */
146  if ( clockwise )
147  {
148  increment *= -1;
149  /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
150  if ( a3 > a1 )
151  a3 -= 2.0 * M_PI;
152  if ( a2 > a1 )
153  a2 -= 2.0 * M_PI;
154  }
155  /* p2 on right side => counter-clockwise sweep */
156  else
157  {
158  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
159  if ( a3 < a1 )
160  a3 += 2.0 * M_PI;
161  if ( a2 < a1 )
162  a2 += 2.0 * M_PI;
163  }
164 
165  /* Override angles for circle case */
166  if( is_circle )
167  {
168  a3 = a1 + 2.0 * M_PI;
169  a2 = a1 + M_PI;
170  increment = fabs(increment);
171  clockwise = LW_FALSE;
172  }
173 
174  /* Initialize point array */
175  pa = ptarray_construct_empty(1, 1, 32);
176 
177  /* Sweep from a1 to a3 */
179  for ( angle = a1 + increment; clockwise ? angle > a3 : angle < a3; angle += increment )
180  {
181  pt.x = center.x + radius * cos(angle);
182  pt.y = center.y + radius * sin(angle);
183  pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
184  pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
185  ptarray_append_point(pa, &pt, LW_FALSE);
186  }
187  return pa;
188 }
189 
190 LWLINE *
191 lwcircstring_segmentize(const LWCIRCSTRING *icurve, uint32_t perQuad)
192 {
193  LWLINE *oline;
194  POINTARRAY *ptarray;
195  POINTARRAY *tmp;
196  uint32_t i, j;
197  POINT4D p1, p2, p3, p4;
198 
199  LWDEBUGF(2, "lwcircstring_segmentize called., dim = %d", icurve->points->flags);
200 
201  ptarray = ptarray_construct_empty(FLAGS_GET_Z(icurve->points->flags), FLAGS_GET_M(icurve->points->flags), 64);
202 
203  for (i = 2; i < icurve->points->npoints; i+=2)
204  {
205  LWDEBUGF(3, "lwcircstring_segmentize: arc ending at point %d", i);
206 
207  getPoint4d_p(icurve->points, i - 2, &p1);
208  getPoint4d_p(icurve->points, i - 1, &p2);
209  getPoint4d_p(icurve->points, i, &p3);
210  tmp = lwcircle_segmentize(&p1, &p2, &p3, perQuad);
211 
212  if (tmp)
213  {
214  LWDEBUGF(3, "lwcircstring_segmentize: generated %d points", tmp->npoints);
215 
216  for (j = 0; j < tmp->npoints; j++)
217  {
218  getPoint4d_p(tmp, j, &p4);
219  ptarray_append_point(ptarray, &p4, LW_TRUE);
220  }
221  ptarray_free(tmp);
222  }
223  else
224  {
225  LWDEBUG(3, "lwcircstring_segmentize: points are colinear, returning curve points as line");
226 
227  for (j = i - 2 ; j < i ; j++)
228  {
229  getPoint4d_p(icurve->points, j, &p4);
230  ptarray_append_point(ptarray, &p4, LW_TRUE);
231  }
232  }
233 
234  }
235  getPoint4d_p(icurve->points, icurve->points->npoints-1, &p1);
236  ptarray_append_point(ptarray, &p1, LW_TRUE);
237 
238  oline = lwline_construct(icurve->srid, NULL, ptarray);
239  return oline;
240 }
241 
242 LWLINE *
243 lwcompound_segmentize(const LWCOMPOUND *icompound, uint32_t perQuad)
244 {
245  LWGEOM *geom;
246  POINTARRAY *ptarray = NULL, *ptarray_out = NULL;
247  LWLINE *tmp = NULL;
248  uint32_t i, j;
249  POINT4D p;
250 
251  LWDEBUG(2, "lwcompound_segmentize called.");
252 
253  ptarray = ptarray_construct_empty(FLAGS_GET_Z(icompound->flags), FLAGS_GET_M(icompound->flags), 64);
254 
255  for (i = 0; i < icompound->ngeoms; i++)
256  {
257  geom = icompound->geoms[i];
258  if (geom->type == CIRCSTRINGTYPE)
259  {
260  tmp = lwcircstring_segmentize((LWCIRCSTRING *)geom, perQuad);
261  for (j = 0; j < tmp->points->npoints; j++)
262  {
263  getPoint4d_p(tmp->points, j, &p);
264  ptarray_append_point(ptarray, &p, LW_TRUE);
265  }
266  lwline_free(tmp);
267  }
268  else if (geom->type == LINETYPE)
269  {
270  tmp = (LWLINE *)geom;
271  for (j = 0; j < tmp->points->npoints; j++)
272  {
273  getPoint4d_p(tmp->points, j, &p);
274  ptarray_append_point(ptarray, &p, LW_TRUE);
275  }
276  }
277  else
278  {
279  lwerror("Unsupported geometry type %d found.",
280  geom->type, lwtype_name(geom->type));
281  return NULL;
282  }
283  }
284  ptarray_out = ptarray_remove_repeated_points(ptarray);
285  ptarray_free(ptarray);
286  return lwline_construct(icompound->srid, NULL, ptarray_out);
287 }
288 
289 LWPOLY *
290 lwcurvepoly_segmentize(const LWCURVEPOLY *curvepoly, uint32_t perQuad)
291 {
292  LWPOLY *ogeom;
293  LWGEOM *tmp;
294  LWLINE *line;
295  POINTARRAY **ptarray;
296  int i;
297 
298  LWDEBUG(2, "lwcurvepoly_segmentize called.");
299 
300  ptarray = lwalloc(sizeof(POINTARRAY *)*curvepoly->nrings);
301 
302  for (i = 0; i < curvepoly->nrings; i++)
303  {
304  tmp = curvepoly->rings[i];
305  if (tmp->type == CIRCSTRINGTYPE)
306  {
307  line = lwcircstring_segmentize((LWCIRCSTRING *)tmp, perQuad);
308  ptarray[i] = ptarray_clone_deep(line->points);
309  lwline_free(line);
310  }
311  else if (tmp->type == LINETYPE)
312  {
313  line = (LWLINE *)tmp;
314  ptarray[i] = ptarray_clone_deep(line->points);
315  }
316  else if (tmp->type == COMPOUNDTYPE)
317  {
318  line = lwcompound_segmentize((LWCOMPOUND *)tmp, perQuad);
319  ptarray[i] = ptarray_clone_deep(line->points);
320  lwline_free(line);
321  }
322  else
323  {
324  lwerror("Invalid ring type found in CurvePoly.");
325  return NULL;
326  }
327  }
328 
329  ogeom = lwpoly_construct(curvepoly->srid, NULL, curvepoly->nrings, ptarray);
330  return ogeom;
331 }
332 
333 LWMLINE *
334 lwmcurve_segmentize(LWMCURVE *mcurve, uint32_t perQuad)
335 {
336  LWMLINE *ogeom;
337  LWGEOM *tmp;
338  LWGEOM **lines;
339  int i;
340 
341  LWDEBUGF(2, "lwmcurve_segmentize called, geoms=%d, dim=%d.", mcurve->ngeoms, FLAGS_NDIMS(mcurve->flags));
342 
343  lines = lwalloc(sizeof(LWGEOM *)*mcurve->ngeoms);
344 
345  for (i = 0; i < mcurve->ngeoms; i++)
346  {
347  tmp = mcurve->geoms[i];
348  if (tmp->type == CIRCSTRINGTYPE)
349  {
350  lines[i] = (LWGEOM *)lwcircstring_segmentize((LWCIRCSTRING *)tmp, perQuad);
351  }
352  else if (tmp->type == LINETYPE)
353  {
354  lines[i] = (LWGEOM *)lwline_construct(mcurve->srid, NULL, ptarray_clone_deep(((LWLINE *)tmp)->points));
355  }
356  else if (tmp->type == COMPOUNDTYPE)
357  {
358  lines[i] = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)tmp, perQuad);
359  }
360  else
361  {
362  lwerror("Unsupported geometry found in MultiCurve.");
363  return NULL;
364  }
365  }
366 
367  ogeom = (LWMLINE *)lwcollection_construct(MULTILINETYPE, mcurve->srid, NULL, mcurve->ngeoms, lines);
368  return ogeom;
369 }
370 
371 LWMPOLY *
372 lwmsurface_segmentize(LWMSURFACE *msurface, uint32_t perQuad)
373 {
374  LWMPOLY *ogeom;
375  LWGEOM *tmp;
376  LWPOLY *poly;
377  LWGEOM **polys;
378  POINTARRAY **ptarray;
379  int i, j;
380 
381  LWDEBUG(2, "lwmsurface_segmentize called.");
382 
383  polys = lwalloc(sizeof(LWGEOM *)*msurface->ngeoms);
384 
385  for (i = 0; i < msurface->ngeoms; i++)
386  {
387  tmp = msurface->geoms[i];
388  if (tmp->type == CURVEPOLYTYPE)
389  {
390  polys[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)tmp, perQuad);
391  }
392  else if (tmp->type == POLYGONTYPE)
393  {
394  poly = (LWPOLY *)tmp;
395  ptarray = lwalloc(sizeof(POINTARRAY *)*poly->nrings);
396  for (j = 0; j < poly->nrings; j++)
397  {
398  ptarray[j] = ptarray_clone_deep(poly->rings[j]);
399  }
400  polys[i] = (LWGEOM *)lwpoly_construct(msurface->srid, NULL, poly->nrings, ptarray);
401  }
402  }
403  ogeom = (LWMPOLY *)lwcollection_construct(MULTIPOLYGONTYPE, msurface->srid, NULL, msurface->ngeoms, polys);
404  return ogeom;
405 }
406 
407 LWCOLLECTION *
408 lwcollection_segmentize(LWCOLLECTION *collection, uint32_t perQuad)
409 {
410  LWCOLLECTION *ocol;
411  LWGEOM *tmp;
412  LWGEOM **geoms;
413  int i;
414 
415  LWDEBUG(2, "lwcollection_segmentize called.");
416 
417  geoms = lwalloc(sizeof(LWGEOM *)*collection->ngeoms);
418 
419  for (i=0; i<collection->ngeoms; i++)
420  {
421  tmp = collection->geoms[i];
422  switch (tmp->type)
423  {
424  case CIRCSTRINGTYPE:
425  geoms[i] = (LWGEOM *)lwcircstring_segmentize((LWCIRCSTRING *)tmp, perQuad);
426  break;
427  case COMPOUNDTYPE:
428  geoms[i] = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)tmp, perQuad);
429  break;
430  case CURVEPOLYTYPE:
431  geoms[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)tmp, perQuad);
432  break;
433  case COLLECTIONTYPE:
434  geoms[i] = (LWGEOM *)lwcollection_segmentize((LWCOLLECTION *)tmp, perQuad);
435  break;
436  default:
437  geoms[i] = lwgeom_clone(tmp);
438  break;
439  }
440  }
441  ocol = lwcollection_construct(COLLECTIONTYPE, collection->srid, NULL, collection->ngeoms, geoms);
442  return ocol;
443 }
444 
445 LWGEOM *
446 lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad)
447 {
448  LWGEOM * ogeom = NULL;
449  switch (geom->type)
450  {
451  case CIRCSTRINGTYPE:
452  ogeom = (LWGEOM *)lwcircstring_segmentize((LWCIRCSTRING *)geom, perQuad);
453  break;
454  case COMPOUNDTYPE:
455  ogeom = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)geom, perQuad);
456  break;
457  case CURVEPOLYTYPE:
458  ogeom = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)geom, perQuad);
459  break;
460  case MULTICURVETYPE:
461  ogeom = (LWGEOM *)lwmcurve_segmentize((LWMCURVE *)geom, perQuad);
462  break;
463  case MULTISURFACETYPE:
464  ogeom = (LWGEOM *)lwmsurface_segmentize((LWMSURFACE *)geom, perQuad);
465  break;
466  case COLLECTIONTYPE:
467  ogeom = (LWGEOM *)lwcollection_segmentize((LWCOLLECTION *)geom, perQuad);
468  break;
469  default:
470  ogeom = lwgeom_clone(geom);
471  }
472  return ogeom;
473 }
474 
479 static double
480 lw_arc_angle(const POINT2D *a, const POINT2D *b, const POINT2D *c)
481 {
482  POINT2D ab, cb;
483 
484  ab.x = b->x - a->x;
485  ab.y = b->y - a->y;
486 
487  cb.x = b->x - c->x;
488  cb.y = b->y - c->y;
489 
490  double dot = (ab.x * cb.x + ab.y * cb.y); /* dot product */
491  double cross = (ab.x * cb.y - ab.y * cb.x); /* cross product */
492 
493  double alpha = atan2(cross, dot);
494 
495  return alpha;
496 }
497 
502 static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b)
503 {
504  POINT2D center;
505  POINT2D *t1 = (POINT2D*)a1;
506  POINT2D *t2 = (POINT2D*)a2;
507  POINT2D *t3 = (POINT2D*)a3;
508  POINT2D *tb = (POINT2D*)b;
509  double radius = lw_arc_center(t1, t2, t3, &center);
510  double b_distance, diff;
511 
512  /* Co-linear a1/a2/a3 */
513  if ( radius < 0.0 )
514  return LW_FALSE;
515 
516  b_distance = distance2d_pt_pt(tb, &center);
517  diff = fabs(radius - b_distance);
518  LWDEBUGF(4, "circle_radius=%g, b_distance=%g, diff=%g, percentage=%g", radius, b_distance, diff, diff/radius);
519 
520  /* Is the point b on the circle? */
521  if ( diff < EPSILON_SQLMM )
522  {
523  int a2_side = lw_segment_side(t1, t3, t2);
524  int b_side = lw_segment_side(t1, t3, tb);
525  double angle1 = lw_arc_angle(t1, t2, t3);
526  double angle2 = lw_arc_angle(t2, t3, tb);
527 
528  /* Is the angle similar to the previous one ? */
529  diff = fabs(angle1 - angle2);
530  LWDEBUGF(4, " angle1: %g, angle2: %g, diff:%g", angle1, angle2, diff);
531  if ( diff > EPSILON_SQLMM )
532  {
533  return LW_FALSE;
534  }
535 
536  /* Is the point b on the same side of a1/a3 as the mid-point a2 is? */
537  /* If not, it's in the unbounded part of the circle, so it continues the arc, return true. */
538  if ( b_side != a2_side )
539  return LW_TRUE;
540  }
541  return LW_FALSE;
542 }
543 
544 static LWGEOM*
545 linestring_from_pa(const POINTARRAY *pa, int srid, int start, int end)
546 {
547  int i = 0, j = 0;
548  POINT4D p;
549  POINTARRAY *pao = ptarray_construct(ptarray_has_z(pa), ptarray_has_m(pa), end-start+2);
550  LWDEBUGF(4, "srid=%d, start=%d, end=%d", srid, start, end);
551  for( i = start; i < end + 2; i++ )
552  {
553  getPoint4d_p(pa, i, &p);
554  ptarray_set_point4d(pao, j++, &p);
555  }
556  return lwline_as_lwgeom(lwline_construct(srid, NULL, pao));
557 }
558 
559 static LWGEOM*
560 circstring_from_pa(const POINTARRAY *pa, int srid, int start, int end)
561 {
562 
563  POINT4D p0, p1, p2;
565  LWDEBUGF(4, "srid=%d, start=%d, end=%d", srid, start, end);
566  getPoint4d_p(pa, start, &p0);
567  ptarray_set_point4d(pao, 0, &p0);
568  getPoint4d_p(pa, (start+end+1)/2, &p1);
569  ptarray_set_point4d(pao, 1, &p1);
570  getPoint4d_p(pa, end+1, &p2);
571  ptarray_set_point4d(pao, 2, &p2);
572  return lwcircstring_as_lwgeom(lwcircstring_construct(srid, NULL, pao));
573 }
574 
575 static LWGEOM*
576 geom_from_pa(const POINTARRAY *pa, int srid, int is_arc, int start, int end)
577 {
578  LWDEBUGF(4, "srid=%d, is_arc=%d, start=%d, end=%d", srid, is_arc, start, end);
579  if ( is_arc )
580  return circstring_from_pa(pa, srid, start, end);
581  else
582  return linestring_from_pa(pa, srid, start, end);
583 }
584 
585 LWGEOM*
586 pta_desegmentize(POINTARRAY *points, int type, int srid)
587 {
588  int i = 0, j, k;
589  POINT4D a1, a2, a3, b;
590  POINT4D first, center;
591  char *edges_in_arcs;
592  int found_arc = LW_FALSE;
593  int current_arc = 1;
594  int num_edges;
595  int edge_type; /* non-zero if edge is part of an arc */
596  int start, end;
597  LWCOLLECTION *outcol;
598  /* Minimum number of edges, per quadrant, required to define an arc */
599  const unsigned int min_quad_edges = 2;
600 
601  /* Die on null input */
602  if ( ! points )
603  lwerror("pta_desegmentize called with null pointarray");
604 
605  /* Null on empty input? */
606  if ( points->npoints == 0 )
607  return NULL;
608 
609  /* We can't desegmentize anything shorter than four points */
610  if ( points->npoints < 4 )
611  {
612  /* Return a linestring here*/
613  lwerror("pta_desegmentize needs implementation for npoints < 4");
614  }
615 
616  /* Allocate our result array of vertices that are part of arcs */
617  num_edges = points->npoints - 1;
618  edges_in_arcs = lwalloc(num_edges + 1);
619  memset(edges_in_arcs, 0, num_edges + 1);
620 
621  /* We make a candidate arc of the first two edges, */
622  /* And then see if the next edge follows it */
623  while( i < num_edges-2 )
624  {
625  unsigned int arc_edges;
626  double num_quadrants;
627  double angle;
628 
629  found_arc = LW_FALSE;
630  /* Make candidate arc */
631  getPoint4d_p(points, i , &a1);
632  getPoint4d_p(points, i+1, &a2);
633  getPoint4d_p(points, i+2, &a3);
634  memcpy(&first, &a1, sizeof(POINT4D));
635 
636  for( j = i+3; j < num_edges+1; j++ )
637  {
638  LWDEBUGF(4, "i=%d, j=%d", i, j);
639  getPoint4d_p(points, j, &b);
640  /* Does this point fall on our candidate arc? */
641  if ( pt_continues_arc(&a1, &a2, &a3, &b) )
642  {
643  /* Yes. Mark this edge and the two preceding it as arc components */
644  LWDEBUGF(4, "pt_continues_arc #%d", current_arc);
645  found_arc = LW_TRUE;
646  for ( k = j-1; k > j-4; k-- )
647  edges_in_arcs[k] = current_arc;
648  }
649  else
650  {
651  /* No. So we're done with this candidate arc */
652  LWDEBUG(4, "pt_continues_arc = false");
653  current_arc++;
654  break;
655  }
656 
657  memcpy(&a1, &a2, sizeof(POINT4D));
658  memcpy(&a2, &a3, sizeof(POINT4D));
659  memcpy(&a3, &b, sizeof(POINT4D));
660  }
661  /* Jump past all the edges that were added to the arc */
662  if ( found_arc )
663  {
664  /* Check if an arc was composed by enough edges to be
665  * really considered an arc
666  * See http://trac.osgeo.org/postgis/ticket/2420
667  */
668  arc_edges = j - 1 - i;
669  LWDEBUGF(4, "arc defined by %d edges found", arc_edges);
670  if ( first.x == b.x && first.y == b.y ) {
671  LWDEBUG(4, "arc is a circle");
672  num_quadrants = 4;
673  }
674  else {
675  lw_arc_center((POINT2D*)&first, (POINT2D*)&b, (POINT2D*)&a1, (POINT2D*)&center);
676  angle = lw_arc_angle((POINT2D*)&first, (POINT2D*)&center, (POINT2D*)&b);
677  int p2_side = lw_segment_side((POINT2D*)&first, (POINT2D*)&a1, (POINT2D*)&b);
678  if ( p2_side >= 0 ) angle = -angle;
679 
680  if ( angle < 0 ) angle = 2 * M_PI + angle;
681  num_quadrants = ( 4 * angle ) / ( 2 * M_PI );
682  LWDEBUGF(4, "arc angle (%g %g, %g %g, %g %g) is %g (side is %d), quandrants:%g", first.x, first.y, center.x, center.y, b.x, b.y, angle, p2_side, num_quadrants);
683  }
684  /* a1 is first point, b is last point */
685  if ( arc_edges < min_quad_edges * num_quadrants ) {
686  LWDEBUGF(4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants);
687  for ( k = j-1; k >= i; k-- )
688  edges_in_arcs[k] = 0;
689  }
690  i = j-1;
691  }
692  else
693  {
694  /* Mark this edge as a linear edge */
695  edges_in_arcs[i] = 0;
696  i = i+1;
697  }
698  }
699 
700 #if POSTGIS_DEBUG_LEVEL > 3
701  {
702  char *edgestr = lwalloc(num_edges+1);
703  for ( i = 0; i < num_edges; i++ )
704  {
705  if ( edges_in_arcs[i] )
706  edgestr[i] = 48 + edges_in_arcs[i];
707  else
708  edgestr[i] = '.';
709  }
710  edgestr[num_edges] = 0;
711  LWDEBUGF(3, "edge pattern %s", edgestr);
712  lwfree(edgestr);
713  }
714 #endif
715 
716  start = 0;
717  edge_type = edges_in_arcs[0];
718  outcol = lwcollection_construct_empty(COMPOUNDTYPE, srid, ptarray_has_z(points), ptarray_has_m(points));
719  for( i = 1; i < num_edges; i++ )
720  {
721  if( edge_type != edges_in_arcs[i] )
722  {
723  end = i - 1;
724  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
725  start = i;
726  edge_type = edges_in_arcs[i];
727  }
728  }
729  lwfree(edges_in_arcs); /* not needed anymore */
730 
731  /* Roll out last item */
732  end = num_edges - 1;
733  lwcollection_add_lwgeom(outcol, geom_from_pa(points, srid, edge_type, start, end));
734 
735  /* Strip down to singleton if only one entry */
736  if ( outcol->ngeoms == 1 )
737  {
738  LWGEOM *outgeom = outcol->geoms[0];
739  outcol->ngeoms = 0; lwcollection_free(outcol);
740  return outgeom;
741  }
742  return lwcollection_as_lwgeom(outcol);
743 }
744 
745 
746 LWGEOM *
748 {
749  LWDEBUG(2, "lwline_desegmentize called.");
750 
751  if ( line->points->npoints < 4 ) return lwline_as_lwgeom(lwline_clone(line));
752  else return pta_desegmentize(line->points, line->flags, line->srid);
753 }
754 
755 LWGEOM *
757 {
758  LWGEOM **geoms;
759  int i, hascurve = 0;
760 
761  LWDEBUG(2, "lwpolygon_desegmentize called.");
762 
763  geoms = lwalloc(sizeof(LWGEOM *)*poly->nrings);
764  for (i=0; i<poly->nrings; i++)
765  {
766  geoms[i] = pta_desegmentize(poly->rings[i], poly->flags, poly->srid);
767  if (geoms[i]->type == CIRCSTRINGTYPE || geoms[i]->type == COMPOUNDTYPE)
768  {
769  hascurve = 1;
770  }
771  }
772  if (hascurve == 0)
773  {
774  for (i=0; i<poly->nrings; i++)
775  {
776  lwfree(geoms[i]); /* TODO: should this be lwgeom_free instead ? */
777  }
778  return lwgeom_clone((LWGEOM *)poly);
779  }
780 
781  return (LWGEOM *)lwcollection_construct(CURVEPOLYTYPE, poly->srid, NULL, poly->nrings, geoms);
782 }
783 
784 LWGEOM *
786 {
787  LWGEOM **geoms;
788  int i, hascurve = 0;
789 
790  LWDEBUG(2, "lwmline_desegmentize called.");
791 
792  geoms = lwalloc(sizeof(LWGEOM *)*mline->ngeoms);
793  for (i=0; i<mline->ngeoms; i++)
794  {
795  geoms[i] = lwline_desegmentize((LWLINE *)mline->geoms[i]);
796  if (geoms[i]->type == CIRCSTRINGTYPE || geoms[i]->type == COMPOUNDTYPE)
797  {
798  hascurve = 1;
799  }
800  }
801  if (hascurve == 0)
802  {
803  for (i=0; i<mline->ngeoms; i++)
804  {
805  lwfree(geoms[i]); /* TODO: should this be lwgeom_free instead ? */
806  }
807  return lwgeom_clone((LWGEOM *)mline);
808  }
809  return (LWGEOM *)lwcollection_construct(MULTICURVETYPE, mline->srid, NULL, mline->ngeoms, geoms);
810 }
811 
812 LWGEOM *
814 {
815  LWGEOM **geoms;
816  int i, hascurve = 0;
817 
818  LWDEBUG(2, "lwmpoly_desegmentize called.");
819 
820  geoms = lwalloc(sizeof(LWGEOM *)*mpoly->ngeoms);
821  for (i=0; i<mpoly->ngeoms; i++)
822  {
823  geoms[i] = lwpolygon_desegmentize((LWPOLY *)mpoly->geoms[i]);
824  if (geoms[i]->type == CURVEPOLYTYPE)
825  {
826  hascurve = 1;
827  }
828  }
829  if (hascurve == 0)
830  {
831  for (i=0; i<mpoly->ngeoms; i++)
832  {
833  lwfree(geoms[i]); /* TODO: should this be lwgeom_free instead ? */
834  }
835  return lwgeom_clone((LWGEOM *)mpoly);
836  }
837  return (LWGEOM *)lwcollection_construct(MULTISURFACETYPE, mpoly->srid, NULL, mpoly->ngeoms, geoms);
838 }
839 
840 LWGEOM *
842 {
843  LWDEBUG(2, "lwgeom_desegmentize called.");
844 
845  switch (geom->type)
846  {
847  case LINETYPE:
848  return lwline_desegmentize((LWLINE *)geom);
849  case POLYGONTYPE:
850  return lwpolygon_desegmentize((LWPOLY *)geom);
851  case MULTILINETYPE:
852  return lwmline_desegmentize((LWMLINE *)geom);
853  case MULTIPOLYGONTYPE:
854  return lwmpolygon_desegmentize((LWMPOLY *)geom);
855  default:
856  return lwgeom_clone(geom);
857  }
858 }
859 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:501
int ngeoms
Definition: liblwgeom.h:502
int32_t srid
Definition: liblwgeom.h:475
double x
Definition: liblwgeom.h:308
double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
Determines the center of the circle defined by the three given points.
Definition: lwalgorithm.c:228
#define LINETYPE
Definition: liblwgeom.h:61
LWPOLY * lwcurvepoly_segmentize(const LWCURVEPOLY *curvepoly, uint32_t perQuad)
Definition: lwsegmentize.c:290
uint8_t flags
Definition: liblwgeom.h:499
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:49
double m
Definition: liblwgeom.h:308
#define MULTICURVETYPE
Definition: liblwgeom.h:70
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
static double lw_arc_angle(const POINT2D *a, const POINT2D *b, const POINT2D *c)
Return ABC angle in radians TODO: move to lwalgorithm.
Definition: lwsegmentize.c:480
void lwfree(void *mem)
Definition: lwutil.c:190
LWGEOM ** rings
Definition: liblwgeom.h:491
LWLINE * lwline_clone(const LWLINE *lwgeom)
Definition: lwline.c:89
int npoints
Definition: liblwgeom.h:327
static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b)
Returns LW_TRUE if b is on the arc formed by a1/a2/a3, but not within that portion already described ...
Definition: lwsegmentize.c:502
int32_t srid
Definition: liblwgeom.h:399
LWGEOM * lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad)
Definition: lwsegmentize.c:446
LWMLINE * lwmcurve_segmentize(LWMCURVE *mcurve, uint32_t perQuad)
Definition: lwsegmentize.c:334
#define POLYGONTYPE
Definition: liblwgeom.h:62
LWGEOM * lwcircstring_as_lwgeom(const LWCIRCSTRING *obj)
Definition: lwgeom.c:224
#define CURVEPOLYTYPE
Definition: liblwgeom.h:69
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:57
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:315
#define COMPOUNDTYPE
Definition: liblwgeom.h:68
#define MULTIPOINTTYPE
Definition: liblwgeom.h:63
void lwline_free(LWLINE *line)
Definition: lwline.c:63
int32_t srid
Definition: liblwgeom.h:501
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
#define TRIANGLETYPE
Definition: liblwgeom.h:73
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:72
static LWGEOM * circstring_from_pa(const POINTARRAY *pa, int srid, int start, int end)
Definition: lwsegmentize.c:560
LWGEOM * lwgeom_desegmentize(LWGEOM *geom)
Definition: lwsegmentize.c:841
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2123
int32_t srid
Definition: liblwgeom.h:377
LWGEOM ** geoms
Definition: liblwgeom.h:478
uint8_t flags
Definition: liblwgeom.h:473
LWGEOM * pta_desegmentize(POINTARRAY *points, int type, int srid)
Definition: lwsegmentize.c:586
int ngeoms
Definition: liblwgeom.h:437
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double x
Definition: liblwgeom.h:284
LWGEOM * lwpolygon_desegmentize(LWPOLY *poly)
Definition: lwsegmentize.c:756
LWLINE * lwcircstring_segmentize(const LWCIRCSTRING *icurve, uint32_t perQuad)
Definition: lwsegmentize.c:191
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
static double interpolate_arc(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
Definition: lwsegmentize.c:81
LWGEOM * lwmline_desegmentize(LWMLINE *mline)
Definition: lwsegmentize.c:785
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_TRUE, then a duplicate point will not be added.
Definition: ptarray.c:141
int32_t srid
Definition: liblwgeom.h:488
#define LW_FALSE
Definition: liblwgeom.h:52
uint8_t flags
Definition: liblwgeom.h:325
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:29
LWPOLY ** geoms
Definition: liblwgeom.h:452
#define EPSILON_SQLMM
Tolerance used to determine equality.
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
LWLINE * lwcompound_segmentize(const LWCOMPOUND *icompound, uint32_t perQuad)
Definition: lwsegmentize.c:243
LWGEOM ** geoms
Definition: liblwgeom.h:465
#define TINTYPE
Definition: liblwgeom.h:74
POINTARRAY ** rings
Definition: liblwgeom.h:413
POINTARRAY * ptarray_clone_deep(const POINTARRAY *ptarray)
Deep clone a pointarray (also clones serialized pointlist)
Definition: ptarray.c:619
static LWGEOM * geom_from_pa(const POINTARRAY *pa, int srid, int is_arc, int start, int end)
Definition: lwsegmentize.c:576
int nrings
Definition: liblwgeom.h:411
int32_t srid
Definition: liblwgeom.h:462
double y
Definition: liblwgeom.h:284
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
double z
Definition: liblwgeom.h:308
LWGEOM * lwgeom_clone(const LWGEOM *lwgeom)
Clone LWGEOM object.
Definition: lwgeom.c:351
int ngeoms
Definition: liblwgeom.h:450
LWCIRCSTRING * lwcircstring_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwcircstring.c:38
LWGEOM * lwline_desegmentize(LWLINE *line)
Definition: lwsegmentize.c:747
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
LWLINE ** geoms
Definition: liblwgeom.h:439
LWGEOM * lwmpolygon_desegmentize(LWMPOLY *mpoly)
Definition: lwsegmentize.c:813
int ngeoms
Definition: liblwgeom.h:515
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:30
#define MULTISURFACETYPE
Definition: liblwgeom.h:71
int32_t srid
Definition: liblwgeom.h:410
LWGEOM ** geoms
Definition: liblwgeom.h:517
int ngeoms
Definition: liblwgeom.h:476
LWGEOM ** geoms
Definition: liblwgeom.h:504
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:107
int32_t srid
Definition: liblwgeom.h:514
uint8_t type
Definition: liblwgeom.h:352
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:316
POINTARRAY * points
Definition: liblwgeom.h:400
uint8_t flags
Definition: liblwgeom.h:408
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:67
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:23
int lwgeom_has_arc(const LWGEOM *geom)
Definition: lwsegmentize.c:42
void * lwalloc(size_t size)
Definition: lwutil.c:175
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:62
double y
Definition: liblwgeom.h:308
POINTARRAY * ptarray_remove_repeated_points(POINTARRAY *in)
Definition: ptarray.c:1401
#define MULTILINETYPE
Definition: liblwgeom.h:64
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
uint8_t flags
Definition: liblwgeom.h:375
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:118
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:174
int32_t srid
Definition: liblwgeom.h:449
static POINTARRAY * lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad)
Definition: lwsegmentize.c:103
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:217
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
LWCOLLECTION * lwcollection_segmentize(LWCOLLECTION *collection, uint32_t perQuad)
Definition: lwsegmentize.c:408
int32_t srid
Definition: liblwgeom.h:436
static LWGEOM * linestring_from_pa(const POINTARRAY *pa, int srid, int start, int end)
Definition: lwsegmentize.c:545
POINTARRAY * points
Definition: liblwgeom.h:378
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:219
LWMPOLY * lwmsurface_segmentize(LWMSURFACE *msurface, uint32_t perQuad)
Definition: lwsegmentize.c:372