PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ lwarc_linearize()

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 
)
static

Segmentize an arc.

Does not add the final vertex

Parameters
toPOINTARRAY to append segmentized vertices to
p1first point defining the arc
p2second point defining the arc
p3third point defining the arc
toltolerance, semantic driven by tolerance_type
tolerance_typesee LW_LINEARIZE_TOLERANCE_TYPE
flagsLW_LINEARIZE_FLAGS
Returns
number of points appended (0 if collinear), or -1 on error (lwerror would be called).

Definition at line 251 of file lwstroke.c.

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, &center);
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  {
318  break;
320  increment = angle_increment_using_max_deviation(tol, radius);
321  break;
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 
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);
499  }
500 
501  if ( reverse )
502  {
503  int i;
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 }
char * s
Definition: cu_in_wkt.c:23
int ptarray_remove_point(POINTARRAY *pa, uint32_t where)
Remove a point from an existing POINTARRAY.
Definition: ptarray.c:251
@ LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE
Tolerance expresses the maximum angle between the radii generating approximation line vertices,...
Definition: liblwgeom.h:2363
@ LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD
Tolerance expresses the number of segments to use for each quarter of circle (quadrant).
Definition: liblwgeom.h:2348
@ LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION
Tolerance expresses the maximum distance between an arbitrary point on the curve and the closest poin...
Definition: liblwgeom.h:2355
#define LW_FALSE
Definition: liblwgeom.h:94
@ LW_LINEARIZE_FLAG_SYMMETRIC
Symmetric linearization means that the output vertices would be the same no matter the order of the p...
Definition: liblwgeom.h:2372
@ LW_LINEARIZE_FLAG_RETAIN_ANGLE
Retain angle instructs the engine to try its best to retain the requested angle between generating ra...
Definition: liblwgeom.h:2392
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:59
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:125
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:319
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,...
Definition: ptarray.c:147
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:93
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:226
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:37
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:62
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:44
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static double angle_increment_using_max_deviation(double max_deviation, double radius)
Definition: lwstroke.c:155
static double interpolate_arc(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
Definition: lwstroke.c:108
static double angle_increment_using_segments_per_quad(double tol)
Definition: lwstroke.c:131
static double angle_increment_using_max_angle(double tol)
Definition: lwstroke.c:221
double y
Definition: liblwgeom.h:390
double x
Definition: liblwgeom.h:390
double m
Definition: liblwgeom.h:414
double x
Definition: liblwgeom.h:414
double z
Definition: liblwgeom.h:414
double y
Definition: liblwgeom.h:414
uint32_t npoints
Definition: liblwgeom.h:427

References angle_increment_using_max_angle(), angle_increment_using_max_deviation(), angle_increment_using_segments_per_quad(), getPoint4d_p(), interpolate_arc(), lw_arc_center(), LW_FALSE, LW_LINEARIZE_FLAG_RETAIN_ANGLE, LW_LINEARIZE_FLAG_SYMMETRIC, LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE, LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION, LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD, lw_segment_side(), LW_TRUE, LWDEBUG, LWDEBUGF, lwerror(), POINT4D::m, POINTARRAY::npoints, ptarray_append_point(), ptarray_construct_empty(), ptarray_free(), ptarray_has_m(), ptarray_has_z(), ptarray_remove_point(), s, POINT2D::x, POINT4D::x, POINT2D::y, POINT4D::y, and POINT4D::z.

Referenced by lwcircstring_linearize().

Here is the call graph for this function:
Here is the caller graph for this function: