PostGIS  2.5.1dev-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 135 of file lwstroke.c.

References 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(), POINT2D::x, POINT4D::x, POINT2D::y, POINT4D::y, and POINT4D::z.

Referenced by lwcircstring_linearize().

139 {
140  POINT2D center;
141  POINT2D *t1 = (POINT2D*)p1;
142  POINT2D *t2 = (POINT2D*)p2;
143  POINT2D *t3 = (POINT2D*)p3;
144  POINT4D pt;
145  int p2_side = 0;
146  int clockwise = LW_TRUE;
147  double radius; /* Arc radius */
148  double increment; /* Angle per segment */
149  double angle_shift = 0;
150  double a1, a2, a3, angle;
151  POINTARRAY *pa = to;
152  int is_circle = LW_FALSE;
153  int points_added = 0;
154  int reverse = 0;
155 
156  LWDEBUG(2, "lwarc_linearize called.");
157 
158  LWDEBUGF(2, " curve is CIRCULARSTRING(%.15g %.15f, %.15f %.15f, %.15f %15f)",
159  t1->x, t1->y, t2->x, t2->y, t3->x, t3->y);
160 
161  p2_side = lw_segment_side(t1, t3, t2);
162 
163  LWDEBUGF(2, " p2 side is %d", p2_side);
164 
165  /* Force counterclockwise scan if SYMMETRIC operation is requsested */
166  if ( p2_side == -1 && flags & LW_LINEARIZE_FLAG_SYMMETRIC )
167  {
168  /* swap p1-p3 */
169  t1 = (POINT2D*)p3;
170  t3 = (POINT2D*)p1;
171  p1 = (POINT4D*)t1;
172  p3 = (POINT4D*)t3;
173  p2_side = 1;
174  reverse = 1;
175  }
176 
177  radius = lw_arc_center(t1, t2, t3, &center);
178  LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
179 
180  /* Matched start/end points imply circle */
181  if ( p1->x == p3->x && p1->y == p3->y )
182  is_circle = LW_TRUE;
183 
184  /* Negative radius signals straight line, p1/p2/p3 are colinear */
185  if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
186  return 0;
187 
188  /* The side of the p1/p3 line that p2 falls on dictates the sweep
189  direction from p1 to p3. */
190  if ( p2_side == -1 )
191  clockwise = LW_TRUE;
192  else
193  clockwise = LW_FALSE;
194 
195  if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD )
196  {{
197  int perQuad = rint(tol);
198  // error out if tol != perQuad ? (not-round)
199  if ( perQuad != tol )
200  {
201  lwerror("lwarc_linearize: segments per quadrant must be an integer value, got %.15g", tol, perQuad);
202  return -1;
203  }
204  if ( perQuad < 1 )
205  {
206  lwerror("lwarc_linearize: segments per quadrant must be at least 1, got %d", perQuad);
207  return -1;
208  }
209  increment = fabs(M_PI_2 / perQuad);
210  LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g (%g degrees)", perQuad, increment, increment*180/M_PI);
211 
212  }}
213  else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION )
214  {{
215  double halfAngle, maxErr;
216  if ( tol <= 0 )
217  {
218  lwerror("lwarc_linearize: max deviation must be bigger than 0, got %.15g", tol);
219  return -1;
220  }
221 
222  /*
223  * Ref: https://en.wikipedia.org/wiki/Sagitta_(geometry)
224  *
225  * An arc "sagitta" (distance between middle point of arc and
226  * middle point of corresponding chord) is defined as:
227  *
228  * sagitta = radius * ( 1 - cos( angle ) );
229  *
230  * We want our sagitta to be at most "tolerance" long,
231  * and we want to find out angle, so we use the inverse
232  * formula:
233  *
234  * tol = radius * ( 1 - cos( angle ) );
235  * 1 - cos( angle ) = tol/radius
236  * - cos( angle ) = tol/radius - 1
237  * cos( angle ) = - tol/radius + 1
238  * angle = acos( 1 - tol/radius )
239  *
240  * Constraints: 1.0 - tol/radius must be between -1 and 1
241  * which means tol must be between 0 and 2 times
242  * the radius, which makes sense as you cannot have a
243  * sagitta bigger than twice the radius!
244  *
245  */
246  maxErr = tol;
247  if ( maxErr > radius * 2 )
248  {
249  maxErr = radius * 2;
250  LWDEBUGF(2, "lwarc_linearize: tolerance %g is too big, "
251  "using arc-max 2 * radius == %g", tol, maxErr);
252  }
253  do {
254  halfAngle = acos( 1.0 - maxErr / radius );
255  /* TODO: avoid a loop here, going rather straight to
256  * a minimum angle value */
257  if ( halfAngle != 0 ) break;
258  LWDEBUGF(2, "lwarc_linearize: tolerance %g is too small for this arc"
259  " to compute approximation angle, doubling it", maxErr);
260  maxErr *= 2;
261  } while(1);
262  increment = 2 * halfAngle;
263  LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)", tol, radius, halfAngle, increment, increment*180/M_PI);
264  }}
265  else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE )
266  {
267  increment = tol;
268  if ( increment <= 0 )
269  {
270  lwerror("lwarc_linearize: max angle must be bigger than 0, got %.15g", tol);
271  return -1;
272  }
273  }
274  else
275  {
276  lwerror("lwarc_linearize: unsupported tolerance type %d", tolerance_type);
277  return LW_FALSE;
278  }
279 
280  /* Angles of each point that defines the arc section */
281  a1 = atan2(p1->y - center.y, p1->x - center.x);
282  a2 = atan2(p2->y - center.y, p2->x - center.x);
283  a3 = atan2(p3->y - center.y, p3->x - center.x);
284 
285  LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
286  a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);
287 
288  if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
289  {{
290  /* Calculate total arc angle, in radians */
291  double angle = clockwise ? a1 - a3 : a3 - a1;
292  if ( angle < 0 ) angle += M_PI * 2;
293  LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg",
294  angle * 180 / M_PI);
295  if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
296  {{
297  /* Number of steps */
298  int steps = trunc(angle / increment);
299  /* Angle reminder */
300  double angle_reminder = angle - ( increment * steps );
301  angle_shift = angle_reminder / 2.0;
302 
303  LWDEBUGF(2, "lwarc_linearize RETAIN_ANGLE operation requested - "
304  "total angle %g, steps %d, increment %g, reminder %g",
305  angle * 180 / M_PI, steps, increment * 180 / M_PI,
306  angle_reminder * 180 / M_PI);
307  }}
308  else
309  {{
310  /* Number of segments in output */
311  int segs = ceil(angle / increment);
312  /* Tweak increment to be regular for all the arc */
313  increment = angle/segs;
314 
315  LWDEBUGF(2, "lwarc_linearize SYMMETRIC operation requested - "
316  "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d - I:%g",
317  angle*180/M_PI, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
318  segs, increment*180/M_PI);
319  }}
320  }}
321 
322  /* p2 on left side => clockwise sweep */
323  if ( clockwise )
324  {
325  LWDEBUG(2, " Clockwise sweep");
326  increment *= -1;
327  angle_shift *= -1;
328  /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
329  if ( a3 > a1 )
330  a3 -= 2.0 * M_PI;
331  if ( a2 > a1 )
332  a2 -= 2.0 * M_PI;
333  }
334  /* p2 on right side => counter-clockwise sweep */
335  else
336  {
337  LWDEBUG(2, " Counterclockwise sweep");
338  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
339  if ( a3 < a1 )
340  a3 += 2.0 * M_PI;
341  if ( a2 < a1 )
342  a2 += 2.0 * M_PI;
343  }
344 
345  /* Override angles for circle case */
346  if( is_circle )
347  {
348  a3 = a1 + 2.0 * M_PI;
349  a2 = a1 + M_PI;
350  increment = fabs(increment);
351  clockwise = LW_FALSE;
352  }
353 
354  LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
355  angle_shift * 180/M_PI, increment * 180/M_PI);
356 
357  if ( reverse ) {{
358  const int capacity = 8; /* TODO: compute exactly ? */
359  pa = ptarray_construct_empty(ptarray_has_z(to), ptarray_has_m(to), capacity);
360  }}
361 
362  /* Sweep from a1 to a3 */
363  if ( ! reverse )
364  {
366  }
367  ++points_added;
368  if ( angle_shift ) angle_shift -= increment;
369  LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
370  a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
371  for ( angle = a1 + increment + angle_shift; clockwise ? angle > a3 : angle < a3; angle += increment )
372  {
373  LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
374  pt.x = center.x + radius * cos(angle);
375  pt.y = center.y + radius * sin(angle);
376  pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
377  pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
378  ptarray_append_point(pa, &pt, LW_FALSE);
379  ++points_added;
380  angle_shift = 0;
381  }
382 
383  if ( reverse ) {{
384  int i;
386  for ( i=pa->npoints; i>0; i-- ) {
387  getPoint4d_p(pa, i-1, &pt);
388  ptarray_append_point(to, &pt, LW_FALSE);
389  }
390  ptarray_free(pa);
391  }}
392 
393  return points_added;
394 }
double x
Definition: liblwgeom.h:354
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
double m
Definition: liblwgeom.h:354
Tolerance expresses the maximum angle between the radii generating approximation line vertices...
Definition: liblwgeom.h:2222
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:328
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
Tolerance expresses the maximum distance between an arbitrary point on the curve and the closest poin...
Definition: liblwgeom.h:2214
double x
Definition: liblwgeom.h:330
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
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
double y
Definition: liblwgeom.h:330
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:354
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
Symmetric linearization means that the output vertices would be the same no matter the order of the p...
Definition: liblwgeom.h:2231
Tolerance expresses the number of segments to use for each quarter of circle (quadrant).
Definition: liblwgeom.h:2207
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
Retain angle instructs the engine to try its best to retain the requested angle between generating ra...
Definition: liblwgeom.h:2251
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
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:354
#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
uint32_t npoints
Definition: liblwgeom.h:373
Here is the call graph for this function:
Here is the caller graph for this function: