PostGIS  2.5.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 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  p2_side = lw_segment_side(t1, t3, t2);
159 
160  /* Force counterclockwise scan if SYMMETRIC operation is requsested */
161  if ( p2_side == -1 && flags & LW_LINEARIZE_FLAG_SYMMETRIC )
162  {
163  /* swap p1-p3 */
164  t1 = (POINT2D*)p3;
165  t3 = (POINT2D*)p1;
166  p1 = (POINT4D*)t1;
167  p3 = (POINT4D*)t3;
168  p2_side = 1;
169  reverse = 1;
170  }
171 
172  radius = lw_arc_center(t1, t2, t3, &center);
173  LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
174 
175  /* Matched start/end points imply circle */
176  if ( p1->x == p3->x && p1->y == p3->y )
177  is_circle = LW_TRUE;
178 
179  /* Negative radius signals straight line, p1/p2/p3 are colinear */
180  if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
181  return 0;
182 
183  /* The side of the p1/p3 line that p2 falls on dictates the sweep
184  direction from p1 to p3. */
185  if ( p2_side == -1 )
186  clockwise = LW_TRUE;
187  else
188  clockwise = LW_FALSE;
189 
190  if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD )
191  {{
192  int perQuad = rint(tol);
193  // error out if tol != perQuad ? (not-round)
194  if ( perQuad != tol )
195  {
196  lwerror("lwarc_linearize: segments per quadrant must be an integer value, got %.15g", tol, perQuad);
197  return -1;
198  }
199  if ( perQuad < 1 )
200  {
201  lwerror("lwarc_linearize: segments per quadrant must be at least 1, got %d", perQuad);
202  return -1;
203  }
204  increment = fabs(M_PI_2 / perQuad);
205  LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g (%g degrees)", perQuad, increment, increment*180/M_PI);
206 
207  }}
208  else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION )
209  {{
210  double halfAngle, maxErr;
211  if ( tol <= 0 )
212  {
213  lwerror("lwarc_linearize: max deviation must be bigger than 0, got %.15g", tol);
214  return -1;
215  }
216 
217  /*
218  * Ref: https://en.wikipedia.org/wiki/Sagitta_(geometry)
219  *
220  * An arc "sagitta" (distance between middle point of arc and
221  * middle point of corresponding chord) is defined as:
222  *
223  * sagitta = radius * ( 1 - cos( angle ) );
224  *
225  * We want our sagitta to be at most "tolerance" long,
226  * and we want to find out angle, so we use the inverse
227  * formula:
228  *
229  * tol = radius * ( 1 - cos( angle ) );
230  * 1 - cos( angle ) = tol/radius
231  * - cos( angle ) = tol/radius - 1
232  * cos( angle ) = - tol/radius + 1
233  * angle = acos( 1 - tol/radius )
234  *
235  * Constraints: 1.0 - tol/radius must be between -1 and 1
236  * which means tol/radius must be between 0 and 2 times
237  * the radius, which makes sense as you cannot have a
238  * sagitta bigger than twice the radius!
239  *
240  */
241  maxErr = tol;
242  if ( maxErr > radius * 2 )
243  {
244  maxErr = radius * 2;
245  LWDEBUGF(2, "lwarc_linearize: tolerance %g is too big, "
246  "using arc-max 2 * radius == %g", tol, maxErr);
247  }
248  halfAngle = acos( 1.0 - maxErr / radius );
249  increment = 2 * halfAngle;
250  LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)", tol, radius, halfAngle, increment, increment*180/M_PI);
251  }}
252  else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE )
253  {
254  increment = tol;
255  if ( increment <= 0 )
256  {
257  lwerror("lwarc_linearize: max angle must be bigger than 0, got %.15g", tol);
258  return -1;
259  }
260  }
261  else
262  {
263  lwerror("lwarc_linearize: unsupported tolerance type %d", tolerance_type);
264  return LW_FALSE;
265  }
266 
267  /* Angles of each point that defines the arc section */
268  a1 = atan2(p1->y - center.y, p1->x - center.x);
269  a2 = atan2(p2->y - center.y, p2->x - center.x);
270  a3 = atan2(p3->y - center.y, p3->x - center.x);
271 
272  LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
273  a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);
274 
275  if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
276  {{
277  /* Calculate total arc angle, in radians */
278  double angle = clockwise ? a1 - a3 : a3 - a1;
279  if ( angle < 0 ) angle += M_PI * 2;
280  LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg",
281  angle * 180 / M_PI);
282  if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
283  {{
284  /* Number of steps */
285  int steps = trunc(angle / increment);
286  /* Angle reminder */
287  double angle_reminder = angle - ( increment * steps );
288  angle_shift = angle_reminder / 2.0;
289 
290  LWDEBUGF(2, "lwarc_linearize RETAIN_ANGLE operation requested - "
291  "total angle %g, steps %d, increment %g, reminder %g",
292  angle * 180 / M_PI, steps, increment * 180 / M_PI,
293  angle_reminder * 180 / M_PI);
294  }}
295  else
296  {{
297  /* Number of segments in output */
298  int segs = ceil(angle / increment);
299  /* Tweak increment to be regular for all the arc */
300  increment = angle/segs;
301 
302  LWDEBUGF(2, "lwarc_linearize SYMMETRIC operation requested - "
303  "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d - I:%g",
304  angle*180/M_PI, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
305  segs, increment*180/M_PI);
306  }}
307  }}
308 
309  /* p2 on left side => clockwise sweep */
310  if ( clockwise )
311  {
312  LWDEBUG(2, " Clockwise sweep");
313  increment *= -1;
314  angle_shift *= -1;
315  /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
316  if ( a3 > a1 )
317  a3 -= 2.0 * M_PI;
318  if ( a2 > a1 )
319  a2 -= 2.0 * M_PI;
320  }
321  /* p2 on right side => counter-clockwise sweep */
322  else
323  {
324  LWDEBUG(2, " Counterclockwise sweep");
325  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
326  if ( a3 < a1 )
327  a3 += 2.0 * M_PI;
328  if ( a2 < a1 )
329  a2 += 2.0 * M_PI;
330  }
331 
332  /* Override angles for circle case */
333  if( is_circle )
334  {
335  a3 = a1 + 2.0 * M_PI;
336  a2 = a1 + M_PI;
337  increment = fabs(increment);
338  clockwise = LW_FALSE;
339  }
340 
341  LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
342  angle_shift * 180/M_PI, increment * 180/M_PI);
343 
344  if ( reverse ) {{
345  const int capacity = 8; /* TODO: compute exactly ? */
346  pa = ptarray_construct_empty(ptarray_has_z(to), ptarray_has_m(to), capacity);
347  }}
348 
349  /* Sweep from a1 to a3 */
350  if ( ! reverse )
351  {
353  }
354  ++points_added;
355  if ( angle_shift ) angle_shift -= increment;
356  LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
357  a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
358  for ( angle = a1 + increment + angle_shift; clockwise ? angle > a3 : angle < a3; angle += increment )
359  {
360  LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
361  pt.x = center.x + radius * cos(angle);
362  pt.y = center.y + radius * sin(angle);
363  pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
364  pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
365  ptarray_append_point(pa, &pt, LW_FALSE);
366  ++points_added;
367  angle_shift = 0;
368  }
369 
370  if ( reverse ) {{
371  int i;
373  for ( i=pa->npoints; i>0; i-- ) {
374  getPoint4d_p(pa, i-1, &pt);
375  ptarray_append_point(to, &pt, LW_FALSE);
376  }
377  ptarray_free(pa);
378  }}
379 
380  return points_added;
381 }
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:227
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: