Logo Search packages:      
Sourcecode: inkscape version File versions

static double Geom::NewtonRaphsonRootFind ( BezierCurve const   Q,
Point const &  P,
double const  u 
) [static]

Use Newton-Raphson iteration to find better root.

Parameters:
Q Current fitted curve
P Digitized point
u Parameter value for "P"
Returns:
Improved u

Definition at line 575 of file bezier-utils.cpp.

References bezier_pt().

Referenced by reparameterize().

{
    assert( 0.0 <= u );
    assert( u <= 1.0 );

    /* Generate control vertices for Q'. */
    Point Q1[3];
    for (unsigned i = 0; i < 3; i++) {
        Q1[i] = 3.0 * ( Q[i+1] - Q[i] );
    }

    /* Generate control vertices for Q''. */
    Point Q2[2];
    for (unsigned i = 0; i < 2; i++) {
        Q2[i] = 2.0 * ( Q1[i+1] - Q1[i] );
    }

    /* Compute Q(u), Q'(u) and Q''(u). */
    Point const Q_u  = bezier_pt(3, Q, u);
    Point const Q1_u = bezier_pt(2, Q1, u);
    Point const Q2_u = bezier_pt(1, Q2, u);

    /* Compute f(u)/f'(u), where f is the derivative wrt u of distsq(u) = 0.5 * the square of the
       distance from P to Q(u).  Here we're using Newton-Raphson to find a stationary point in the
       distsq(u), hopefully corresponding to a local minimum in distsq (and hence a local minimum
       distance from P to Q(u)). */
    Point const diff = Q_u - P;
    double numerator = dot(diff, Q1_u);
    double denominator = dot(Q1_u, Q1_u) + dot(diff, Q2_u);

    double improved_u;
    if ( denominator > 0. ) {
        /* One iteration of Newton-Raphson:
           improved_u = u - f(u)/f'(u) */
        improved_u = u - ( numerator / denominator );
    } else {
        /* Using Newton-Raphson would move in the wrong direction (towards a local maximum rather
           than local minimum), so we move an arbitrary amount in the right direction. */
        if ( numerator > 0. ) {
            improved_u = u * .98 - .01;
        } else if ( numerator < 0. ) {
            /* Deliberately asymmetrical, to reduce the chance of cycling. */
            improved_u = .031 + u * .98;
        } else {
            improved_u = u;
        }
    }

    if (!IS_FINITE(improved_u)) {
        improved_u = u;
    } else if ( improved_u < 0.0 ) {
        improved_u = 0.0;
    } else if ( improved_u > 1.0 ) {
        improved_u = 1.0;
    }

    /* Ensure that improved_u isn't actually worse. */
    {
        double const diff_lensq = lensq(diff);
        for (double proportion = .125; ; proportion += .125) {
            if ( lensq( bezier_pt(3, Q, improved_u) - P ) > diff_lensq ) {
                if ( proportion > 1.0 ) {
                    //g_warning("found proportion %g", proportion);
                    improved_u = u;
                    break;
                }
                improved_u = ( ( 1 - proportion ) * improved_u  +
                               proportion         * u            );
            } else {
                break;
            }
        }
    }

    DOUBLE_ASSERT(improved_u);
    return improved_u;
}


Generated by  Doxygen 1.6.0   Back to index