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

to	POINTARRAY to append segmentized vertices to
p1	first point defining the arc
p2	second point defining the arc
p3	third point defining the arc
tol	tolerance, semantic driven by tolerance_type
tolerance_type	see LW_LINEARIZE_TOLERANCE_TYPE
flags	LW_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.

{
        POINT2D center;
        POINT2D *t1 = (POINT2D*)p1;
        POINT2D *t2 = (POINT2D*)p2;
        POINT2D *t3 = (POINT2D*)p3;
        POINT4D pt;
        int p2_side = 0;
        int clockwise = LW_TRUE;
        double radius; /* Arc radius */
        double increment; /* Angle per segment */
        double angle_shift = 0;
        double a1, a2, a3, angle;
        POINTARRAY *pa = to;
        int is_circle = LW_FALSE;
        int points_added = 0;
        int reverse = 0;
 
        LWDEBUG(2, "lwarc_linearize called.");
 
        LWDEBUGF(2, " curve is CIRCULARSTRING(%.15g %.15f, %.15f %.15f, %.15f %15f)",
                t1->x, t1->y, t2->x, t2->y, t3->x, t3->y);
 
        p2_side = lw_segment_side(t1, t3, t2);
 
        LWDEBUGF(2, " p2 side is %d", p2_side);
 
        /* Force counterclockwise scan if SYMMETRIC operation is requsested */
        if ( p2_side == -1 && flags & LW_LINEARIZE_FLAG_SYMMETRIC )
        {
                /* swap p1-p3 */
                t1 = (POINT2D*)p3;
                t3 = (POINT2D*)p1;
                p1 = (POINT4D*)t1;
                p3 = (POINT4D*)t3;
                p2_side = 1;
                reverse = 1;
        }
 
        radius = lw_arc_center(t1, t2, t3, &center);
        LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
 
        /* Matched start/end points imply circle */
        if ( p1->x == p3->x && p1->y == p3->y )
                is_circle = LW_TRUE;
 
        /* Negative radius signals straight line, p1/p2/p3 are colinear */
        if ( (radius < 0.0 || p2_side == 0) && ! is_circle )
            return 0;
 
        /* The side of the p1/p3 line that p2 falls on dictates the sweep
           direction from p1 to p3. */
        if ( p2_side == -1 )
                clockwise = LW_TRUE;
        else
                clockwise = LW_FALSE;
 
        if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_SEGS_PER_QUAD )
        {{
                int perQuad = rint(tol);
                // error out if tol != perQuad ? (not-round)
                if ( perQuad != tol )
                {
                        lwerror("lwarc_linearize: segments per quadrant must be an integer value, got %.15g", tol, perQuad);
                        return -1;
                }
                if ( perQuad < 1 )
                {
                        lwerror("lwarc_linearize: segments per quadrant must be at least 1, got %d", perQuad);
                        return -1;
                }
                increment = fabs(M_PI_2 / perQuad);
                LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g (%g degrees)", perQuad, increment, increment*180/M_PI);
 
        }}
        else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION )
        {{
                double halfAngle, maxErr;
                if ( tol <= 0 )
                {
                        lwerror("lwarc_linearize: max deviation must be bigger than 0, got %.15g", tol);
                        return -1;
                }
 
                /*
                 * Ref: https://en.wikipedia.org/wiki/Sagitta_(geometry)
                 *
                 * An arc "sagitta" (distance between middle point of arc and
                 * middle point of corresponding chord) is defined as:
                 *
                 *   sagitta = radius * ( 1 - cos( angle ) );
                 *
                 * We want our sagitta to be at most "tolerance" long,
                 * and we want to find out angle, so we use the inverse
                 * formula:
                 *
                 *   tol = radius * ( 1 - cos( angle ) );
                 *   1 - cos( angle ) =  tol/radius
                 *   - cos( angle ) =  tol/radius - 1
                 *   cos( angle ) =  - tol/radius + 1
                 *   angle = acos( 1 - tol/radius )
                 *
                 * Constraints: 1.0 - tol/radius must be between -1 and 1
                 * which means tol must be between 0 and 2 times
                 * the radius, which makes sense as you cannot have a
                 * sagitta bigger than twice the radius!
                 *
                 */
                maxErr = tol;
                if ( maxErr > radius * 2 )
                {
                        maxErr = radius * 2;
                        LWDEBUGF(2, "lwarc_linearize: tolerance %g is too big, "
                                    "using arc-max 2 * radius == %g", tol, maxErr);
                }
                do {
                        halfAngle = acos( 1.0 - maxErr / radius );
                        /* TODO: avoid a loop here, going rather straight to
                         *       a minimum angle value */
                        if ( halfAngle != 0 ) break;
                        LWDEBUGF(2, "lwarc_linearize: tolerance %g is too small for this arc"
                                                                        " to compute approximation angle, doubling it", maxErr);
                        maxErr *= 2;
                } while(1);
                increment = 2 * halfAngle;
                LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)", tol, radius, halfAngle, increment, increment*180/M_PI);
        }}
        else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE )
        {
                increment = tol;
                if ( increment <= 0 )
                {
                        lwerror("lwarc_linearize: max angle must be bigger than 0, got %.15g", tol);
                        return -1;
                }
        }
        else
        {
                lwerror("lwarc_linearize: unsupported tolerance type %d", tolerance_type);
                return LW_FALSE;
        }
 
        /* Angles of each point that defines the arc section */
        a1 = atan2(p1->y - center.y, p1->x - center.x);
        a2 = atan2(p2->y - center.y, p2->x - center.x);
        a3 = atan2(p3->y - center.y, p3->x - center.x);
 
        LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
                a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);
 
        if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
        {{
                /* Calculate total arc angle, in radians */
                double angle = clockwise ? a1 - a3 : a3 - a1;
                if ( angle < 0 ) angle += M_PI * 2;
                LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg",
                                 angle * 180 / M_PI);
                if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
                {{
                        /* Number of steps */
                        int steps = trunc(angle / increment);
                        /* Angle reminder */
                        double angle_reminder = angle - ( increment * steps );
                        angle_shift = angle_reminder / 2.0;
 
                        LWDEBUGF(2, "lwarc_linearize RETAIN_ANGLE operation requested - "
                                 "total angle %g, steps %d, increment %g, reminder %g",
                                 angle * 180 / M_PI, steps, increment * 180 / M_PI,
                                 angle_reminder * 180 / M_PI);
                }}
                else
                {{
                        /* Number of segments in output */
                        int segs = ceil(angle / increment);
                        /* Tweak increment to be regular for all the arc */
                        increment = angle/segs;
 
                        LWDEBUGF(2, "lwarc_linearize SYMMETRIC operation requested - "
                                                        "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d -   I:%g",
                                                        angle*180/M_PI, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
                                                        segs, increment*180/M_PI);
                }}
        }}
 
        /* p2 on left side => clockwise sweep */
        if ( clockwise )
        {
                LWDEBUG(2, " Clockwise sweep");
                increment *= -1;
                angle_shift *= -1;
                /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
                if ( a3 > a1 )
                        a3 -= 2.0 * M_PI;
                if ( a2 > a1 )
                        a2 -= 2.0 * M_PI;
        }
        /* p2 on right side => counter-clockwise sweep */
        else
        {
                LWDEBUG(2, " Counterclockwise sweep");
                /* Adjust a3 up so we can increment from a1 to a3 cleanly */
                if ( a3 < a1 )
                        a3 += 2.0 * M_PI;
                if ( a2 < a1 )
                        a2 += 2.0 * M_PI;
        }
 
        /* Override angles for circle case */
        if( is_circle )
        {
                increment = fabs(increment);
                a3 = a1 + 2.0 * M_PI;
                a2 = a1 + M_PI;
                clockwise = LW_FALSE;
        }
 
        LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
                angle_shift * 180/M_PI, increment * 180/M_PI);
 
        if ( reverse ) {{
                const int capacity = 8; /* TODO: compute exactly ? */
                pa = ptarray_construct_empty(ptarray_has_z(to), ptarray_has_m(to), capacity);
        }}
 
        /* Sweep from a1 to a3 */
        if ( ! reverse )
        {
                ptarray_append_point(pa, p1, LW_FALSE);
        }
        ++points_added;
        if ( angle_shift ) angle_shift -= increment;
        LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
                a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
        for ( angle = a1 + increment + angle_shift; clockwise ? angle > a3 : angle < a3; angle += increment )
        {
                LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
                pt.x = center.x + radius * cos(angle);
                pt.y = center.y + radius * sin(angle);
                pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);
                pt.m = interpolate_arc(angle, a1, a2, a3, p1->m, p2->m, p3->m);
                ptarray_append_point(pa, &pt, LW_FALSE);
                ++points_added;
                angle_shift = 0;
        }
 
        /* Ensure the final point is EXACTLY the same as the first for the circular case */
        if ( is_circle )
        {
                ptarray_remove_point(pa, pa->npoints - 1);
                ptarray_append_point(pa, p1, LW_FALSE);
        }
 
        if ( reverse )
        {
                int i;
                ptarray_append_point(to, p3, LW_FALSE);
                for ( i=pa->npoints; i>0; i-- ) {
                        getPoint4d_p(pa, i-1, &pt);
                        ptarray_append_point(to, &pt, LW_FALSE);
                }
                ptarray_free(pa);
        }
 
        return points_added;
}

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(), ptarray_remove_point(), 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: