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