PostGIS  2.5.0dev-r@@SVN_REVISION@@
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 134 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().

138 {
139  POINT2D center;
140  POINT2D *t1 = (POINT2D*)p1;
141  POINT2D *t2 = (POINT2D*)p2;
142  POINT2D *t3 = (POINT2D*)p3;
143  POINT4D pt;
144  int p2_side = 0;
145  int clockwise = LW_TRUE;
146  double radius; /* Arc radius */
147  double increment; /* Angle per segment */
148  double angle_shift = 0;
149  double a1, a2, a3, angle;
150  POINTARRAY *pa = to;
151  int is_circle = LW_FALSE;
152  int points_added = 0;
153  int reverse = 0;
154 
155  LWDEBUG(2, "lwarc_linearize called.");
156 
157  p2_side = lw_segment_side(t1, t3, t2);
158 
159  /* Force counterclockwise scan if SYMMETRIC operation is requsested */
160  if ( p2_side == -1 && flags & LW_LINEARIZE_FLAG_SYMMETRIC )
161  {
162  /* swap p1-p3 */
163  t1 = (POINT2D*)p3;
164  t3 = (POINT2D*)p1;
165  p1 = (POINT4D*)t1;
166  p3 = (POINT4D*)t3;
167  p2_side = 1;
168  reverse = 1;
169  }
170 
171  radius = lw_arc_center(t1, t2, t3, &center);
172  LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
173 
174  /* Matched start/end points imply circle */
175  if ( p1->x == p3->x && p1->y == p3->y )
176  is_circle = LW_TRUE;
177 
178  /* Negative radius signals straight line, p1/p2/p3 are colinear */
179  if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
180  return 0;
181 
182  /* The side of the p1/p3 line that p2 falls on dictates the sweep
183  direction from p1 to p3. */
184  if ( p2_side == -1 )
185  clockwise = LW_TRUE;
186  else
187  clockwise = LW_FALSE;
188 
189  if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD )
190  {{
191  int perQuad = rint(tol);
192  // error out if tol != perQuad ? (not-round)
193  if ( perQuad != tol )
194  {
195  lwerror("lwarc_linearize: segments per quadrant must be an integer value, got %.15g", tol, perQuad);
196  return -1;
197  }
198  if ( perQuad < 1 )
199  {
200  lwerror("lwarc_linearize: segments per quadrant must be at least 1, got %d", perQuad);
201  return -1;
202  }
203  increment = fabs(M_PI_2 / perQuad);
204  LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g (%g degrees)", perQuad, increment, increment*180/M_PI);
205 
206  }}
207  else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION )
208  {{
209  double halfAngle;
210  if ( tol <= 0 )
211  {
212  lwerror("lwarc_linearize: max deviation must be bigger than 0, got %.15g", tol);
213  return -1;
214  }
215  halfAngle = acos( -tol / radius + 1 );
216  increment = 2 * halfAngle;
217  LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)", tol, radius, halfAngle, increment, increment*180/M_PI);
218  }}
219  else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE )
220  {
221  increment = tol;
222  if ( increment <= 0 )
223  {
224  lwerror("lwarc_linearize: max angle must be bigger than 0, got %.15g", tol);
225  return -1;
226  }
227  }
228  else
229  {
230  lwerror("lwarc_linearize: unsupported tolerance type %d", tolerance_type);
231  return LW_FALSE;
232  }
233 
234  /* Angles of each point that defines the arc section */
235  a1 = atan2(p1->y - center.y, p1->x - center.x);
236  a2 = atan2(p2->y - center.y, p2->x - center.x);
237  a3 = atan2(p3->y - center.y, p3->x - center.x);
238 
239  LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
240  a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);
241 
242  if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
243  {{
244  /* Calculate total arc angle, in radians */
245  double angle = clockwise ? a1 - a3 : a3 - a1;
246  if ( angle < 0 ) angle += M_PI * 2;
247  LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg",
248  angle * 180 / M_PI);
249  if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
250  {{
251  /* Number of steps */
252  int steps = trunc(angle / increment);
253  /* Angle reminder */
254  double angle_reminder = angle - ( increment * steps );
255  angle_shift = angle_reminder / 2.0;
256 
257  LWDEBUGF(2, "lwarc_linearize RETAIN_ANGLE operation requested - "
258  "total angle %g, steps %d, increment %g, reminder %g",
259  angle * 180 / M_PI, steps, increment * 180 / M_PI,
260  angle_reminder * 180 / M_PI);
261  }}
262  else
263  {{
264  /* Number of segments in output */
265  int segs = ceil(angle / increment);
266  /* Tweak increment to be regular for all the arc */
267  increment = angle/segs;
268 
269  LWDEBUGF(2, "lwarc_linearize SYMMETRIC operation requested - "
270  "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d - I:%g",
271  angle*180/M_PI, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
272  segs, increment*180/M_PI);
273  }}
274  }}
275 
276  /* p2 on left side => clockwise sweep */
277  if ( clockwise )
278  {
279  LWDEBUG(2, " Clockwise sweep");
280  increment *= -1;
281  angle_shift *= -1;
282  /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
283  if ( a3 > a1 )
284  a3 -= 2.0 * M_PI;
285  if ( a2 > a1 )
286  a2 -= 2.0 * M_PI;
287  }
288  /* p2 on right side => counter-clockwise sweep */
289  else
290  {
291  LWDEBUG(2, " Counterclockwise sweep");
292  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
293  if ( a3 < a1 )
294  a3 += 2.0 * M_PI;
295  if ( a2 < a1 )
296  a2 += 2.0 * M_PI;
297  }
298 
299  /* Override angles for circle case */
300  if( is_circle )
301  {
302  a3 = a1 + 2.0 * M_PI;
303  a2 = a1 + M_PI;
304  increment = fabs(increment);
305  clockwise = LW_FALSE;
306  }
307 
308  LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
309  angle_shift * 180/M_PI, increment * 180/M_PI);
310 
311  if ( reverse ) {{
312  const int capacity = 8; /* TODO: compute exactly ? */
313  pa = ptarray_construct_empty(ptarray_has_z(to), ptarray_has_m(to), capacity);
314  }}
315 
316  /* Sweep from a1 to a3 */
317  if ( ! reverse )
318  {
320  }
321  ++points_added;
322  if ( angle_shift ) angle_shift -= increment;
323  LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
324  a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
325  for ( angle = a1 + increment + angle_shift; clockwise ? angle > a3 : angle < a3; angle += increment )
326  {
327  LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
328  pt.x = center.x + radius * cos(angle);
329  pt.y = center.y + radius * sin(angle);
330  pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
331  pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
332  ptarray_append_point(pa, &pt, LW_FALSE);
333  ++points_added;
334  angle_shift = 0;
335  }
336 
337  if ( reverse ) {{
338  int i;
340  for ( i=pa->npoints; i>0; i-- ) {
341  getPoint4d_p(pa, i-1, &pt);
342  ptarray_append_point(to, &pt, LW_FALSE);
343  }
344  ptarray_free(pa);
345  }}
346 
347  return points_added;
348 }
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
double m
Definition: liblwgeom.h:351
Tolerance expresses the maximum angle between the radii generating approximation line vertices...
Definition: liblwgeom.h:2202
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:2194
double x
Definition: liblwgeom.h:327
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:327
static double interpolate_arc(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
Definition: lwstroke.c:95
double z
Definition: liblwgeom.h:351
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:2211
Tolerance expresses the number of segments to use for each quarter of circle (quadrant).
Definition: liblwgeom.h:2187
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:2231
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:351
#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:370

Here is the call graph for this function:

Here is the caller graph for this function: