PostGIS  3.0.6dev-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 239 of file lwstroke.c.

243 {
244  POINT2D center;
245  POINT2D *t1 = (POINT2D*)p1;
246  POINT2D *t2 = (POINT2D*)p2;
247  POINT2D *t3 = (POINT2D*)p3;
248  POINT4D pt;
249  int p2_side = 0;
250  int clockwise = LW_TRUE;
251  double radius; /* Arc radius */
252  double increment; /* Angle per segment */
253  double angle_shift = 0;
254  double a1, a2, a3;
255  POINTARRAY *pa;
256  int is_circle = LW_FALSE;
257  int points_added = 0;
258  int reverse = 0;
259  int segments = 0;
260 
261  LWDEBUG(2, "lwarc_linearize called.");
262 
263  LWDEBUGF(2, " curve is CIRCULARSTRING(%.15g %.15f, %.15f %.15f, %.15f %15f)",
264  t1->x, t1->y, t2->x, t2->y, t3->x, t3->y);
265 
266  p2_side = lw_segment_side(t1, t3, t2);
267 
268  LWDEBUGF(2, " p2 side is %d", p2_side);
269 
270  /* Force counterclockwise scan if SYMMETRIC operation is requested */
271  if ( p2_side == -1 && flags & LW_LINEARIZE_FLAG_SYMMETRIC )
272  {
273  /* swap p1-p3 */
274  t1 = (POINT2D*)p3;
275  t3 = (POINT2D*)p1;
276  p1 = (POINT4D*)t1;
277  p3 = (POINT4D*)t3;
278  p2_side = 1;
279  reverse = 1;
280  }
281 
282  radius = lw_arc_center(t1, t2, t3, &center);
283  LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
284 
285  /* Matched start/end points imply circle */
286  if ( p1->x == p3->x && p1->y == p3->y )
287  is_circle = LW_TRUE;
288 
289  /* Negative radius signals straight line, p1/p2/p3 are collinear */
290  if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
291  return 0;
292 
293  /* The side of the p1/p3 line that p2 falls on dictates the sweep
294  direction from p1 to p3. */
295  if ( p2_side == -1 )
296  clockwise = LW_TRUE;
297  else
298  clockwise = LW_FALSE;
299 
300  /* Compute the increment (angle per segment) depending on
301  * our tolerance type. */
302  switch(tolerance_type)
303  {
306  break;
308  increment = angle_increment_using_max_deviation(tol, radius);
309  break;
311  increment = angle_increment_using_max_angle(tol);
312  break;
313  default:
314  lwerror("lwarc_linearize: unsupported tolerance type %d", tolerance_type);
315  return -1;
316  }
317 
318  if (increment < 0)
319  {
320  /* Error occurred in increment calculation somewhere
321  * (lwerror already called)
322  */
323  return -1;
324  }
325 
326  /* Angles of each point that defines the arc section */
327  a1 = atan2(p1->y - center.y, p1->x - center.x);
328  a2 = atan2(p2->y - center.y, p2->x - center.x);
329  a3 = atan2(p3->y - center.y, p3->x - center.x);
330 
331  LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
332  a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);
333 
334  /* Calculate total arc angle, in radians */
335  double total_angle = clockwise ? a1 - a3 : a3 - a1;
336  if ( total_angle <= 0 ) total_angle += M_PI * 2;
337 
338  /* At extreme tolerance values (very low or very high, depending on
339  * the semantic) we may cause our arc to collapse. In this case,
340  * we want shrink the increment enough so that we get two segments
341  * for a standard arc, or three segments for a complete circle. */
342  int min_segs = is_circle ? 3 : 2;
343  segments = ceil(total_angle / increment);
344  if (segments < min_segs)
345  {
346  segments = min_segs;
347  increment = total_angle / min_segs;
348  }
349 
350  if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
351  {
352  LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg", total_angle * 180 / M_PI);
353 
354  if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
355  {
356  /* Number of complete steps */
357  segments = trunc(total_angle / increment);
358 
359  /* Figure out the angle remainder, i.e. the amount of the angle
360  * that is left after we can take no more complete angle
361  * increments. */
362  double angle_remainder = total_angle - (increment * segments);
363 
364  /* Shift the starting angle by half of the remainder. This
365  * will have the effect of evenly distributing the remainder
366  * among the first and last segments in the arc. */
367  angle_shift = angle_remainder / 2.0;
368 
369  LWDEBUGF(2,
370  "lwarc_linearize RETAIN_ANGLE operation requested - "
371  "total angle %g, steps %d, increment %g, remainder %g",
372  total_angle * 180 / M_PI,
373  segments,
374  increment * 180 / M_PI,
375  angle_remainder * 180 / M_PI);
376  }
377  else
378  {
379  /* Number of segments in output */
380  segments = ceil(total_angle / increment);
381  /* Tweak increment to be regular for all the arc */
382  increment = total_angle / segments;
383 
384  LWDEBUGF(2,
385  "lwarc_linearize SYMMETRIC operation requested - "
386  "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d - I:%g",
387  total_angle * 180 / M_PI,
388  p1->x,
389  p1->y,
390  center.x,
391  center.y,
392  p3->x,
393  p3->y,
394  segments,
395  increment * 180 / M_PI);
396  }
397  }
398 
399  /* p2 on left side => clockwise sweep */
400  if ( clockwise )
401  {
402  LWDEBUG(2, " Clockwise sweep");
403  increment *= -1;
404  angle_shift *= -1;
405  /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
406  if ( a3 > a1 )
407  a3 -= 2.0 * M_PI;
408  if ( a2 > a1 )
409  a2 -= 2.0 * M_PI;
410  }
411  /* p2 on right side => counter-clockwise sweep */
412  else
413  {
414  LWDEBUG(2, " Counterclockwise sweep");
415  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
416  if ( a3 < a1 )
417  a3 += 2.0 * M_PI;
418  if ( a2 < a1 )
419  a2 += 2.0 * M_PI;
420  }
421 
422  /* Override angles for circle case */
423  if( is_circle )
424  {
425  increment = fabs(increment);
426  segments = ceil(total_angle / increment);
427  if (segments < 3)
428  {
429  segments = 3;
430  increment = total_angle / 3;
431  }
432  a3 = a1 + 2.0 * M_PI;
433  a2 = a1 + M_PI;
434  clockwise = LW_FALSE;
435  angle_shift = 0.0;
436  }
437 
438  LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
439  angle_shift * 180/M_PI, increment * 180/M_PI);
440 
441  if ( reverse )
442  {
443  /* Append points in order to a temporary POINTARRAY and
444  * reverse them before writing to the output POINTARRAY. */
445  const int capacity = 8; /* TODO: compute exactly ? */
446  pa = ptarray_construct_empty(ptarray_has_z(to), ptarray_has_m(to), capacity);
447  }
448  else
449  {
450  /* Append points directly to the output POINTARRAY,
451  * starting with p1. */
452  pa = to;
453 
455  ++points_added;
456  }
457 
458  /* Sweep from a1 to a3 */
459  int seg_start = 1; /* First point is added manually */
460  int seg_end = segments;
461  if (angle_shift != 0.0)
462  {
463  /* When we have extra angles we need to add the extra segments at the
464  * start and end that cover those parts of the arc */
465  seg_start = 0;
466  seg_end = segments + 1;
467  }
468  LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
469  a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
470  for (int s = seg_start; s < seg_end; s++)
471  {
472  double angle = a1 + increment * s + angle_shift;
473  LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
474  pt.x = center.x + radius * cos(angle);
475  pt.y = center.y + radius * sin(angle);
476  pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
477  pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
478  ptarray_append_point(pa, &pt, LW_FALSE);
479  ++points_added;
480  }
481 
482  /* Ensure the final point is EXACTLY the same as the first for the circular case */
483  if ( is_circle )
484  {
485  ptarray_remove_point(pa, pa->npoints - 1);
487  }
488 
489  if ( reverse )
490  {
491  int i;
493  for ( i=pa->npoints; i>0; i-- ) {
494  getPoint4d_p(pa, i-1, &pt);
495  ptarray_append_point(to, &pt, LW_FALSE);
496  }
497  ptarray_free(pa);
498  }
499 
500  return points_added;
501 }
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:2278
@ LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD
Tolerance expresses the number of segments to use for each quarter of circle (quadrant).
Definition: liblwgeom.h:2263
@ 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:2270
#define LW_FALSE
Definition: liblwgeom.h:108
@ 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:2287
@ 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:2307
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:107
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:229
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:65
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:143
static double interpolate_arc(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
Definition: lwstroke.c:96
static double angle_increment_using_segments_per_quad(double tol)
Definition: lwstroke.c:119
static double angle_increment_using_max_angle(double tol)
Definition: lwstroke.c:209
double y
Definition: liblwgeom.h:376
double x
Definition: liblwgeom.h:376
double m
Definition: liblwgeom.h:400
double x
Definition: liblwgeom.h:400
double z
Definition: liblwgeom.h:400
double y
Definition: liblwgeom.h:400
uint32_t npoints
Definition: liblwgeom.h:413

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: