C++ Bezier::at1方法代码示例

double Bernsteins::secant(Bezier bz) {
    double s = 0, t = 1;
    double e = 1e-14;
    int side = 0;
    double r, fr, fs = bz.at0(), ft = bz.at1();

    for (size_t n = 0; n < 100; ++n)
    {
        r = (fs*t - ft*s) / (fs - ft);
        if (fabs(t-s) < e * fabs(t+s)) {
            debug(std::cout << "error small " << fabs(t-s) 
                      << ", accepting solution " << r 
                  << "after " << n << "iterations\n");
            return r;
        }

        fr = horner(bz, r);

        if (fr * ft > 0)
        {
            t = r; ft = fr;
            if (side == -1) fs /= 2;
            side = -1;
        }
        else if (fs * fr > 0)
        {
            s = r;  fs = fr;
            if (side == +1) ft /= 2;
            side = +1;
        }
        else break;
    }
    return r;
}

// suggested by Sederberg.
double Bernsteins::horner(Bezier bz, double t)
{
    double u, tn, tmp;
    u = 1.0 - t;
    tn = 1.0;
    tmp = bz.at0() * u;
    for(size_t i = 1; i < bz.degree(); ++i)
    {
        tn *= t;
        tmp = (tmp + tn*choose<double>(bz.order(), (unsigned)i)*bz[i]) * u;
    }
    return (tmp + tn*t*bz.at1());
}

vector<double> find_all_roots(Bezier b) {
    vector<double> rts = b.roots();
    if(b.at0() == 0) rts.push_back(0);
    if(b.at1() == 0) rts.push_back(1);
    return rts;
}

void Bernsteins::find_bernstein_roots(Bezier bz,
                                      unsigned depth,
                                      double left_t,
                                      double right_t)
{
    debug(std::cout << left_t << ", " << right_t << std::endl);
    size_t n_crossings = 0;

    int old_sign = SGN(bz[0]);
    //std::cout << "w[0] = " << bz[0] << std::endl;
    int sign;
    for (size_t i = 1; i < bz.size(); i++)
    {
        //std::cout << "w[" << i << "] = " << w[i] << std::endl;
        sign = SGN(bz[i]);
        if (sign != 0)
        {
            if (sign != old_sign && old_sign != 0)
            {
                ++n_crossings;
            }
            old_sign = sign;
        }
    }
    //std::cout << "n_crossings = " << n_crossings << std::endl;
    if (n_crossings == 0)  return; // no solutions here

    if (n_crossings == 1) /* Unique solution  */
    {
        //std::cout << "depth = " << depth << std::endl;
        /* Stop recursion when the tree is deep enough  */
        /* if deep enough, return 1 solution at midpoint  */
        if (depth > MAX_DEPTH)
        {
            //printf("bottom out %d\n", depth);
            const double Ax = right_t - left_t;
            const double Ay = bz.at1() - bz.at0();

            solutions.push_back(left_t - Ax*bz.at0() / Ay);
            return;
        }

        double r = secant(bz);
        solutions.push_back(r*right_t + (1-r)*left_t);
        return;
    }
    /* Otherwise, solve recursively after subdividing control polygon  */
    Bezier::Order o(bz);
    Bezier Left(o), Right = bz;
    double split_t = (left_t + right_t) * 0.5;

    // If subdivision is working poorly, split around the leftmost root of the derivative
    if (depth > 2) {
        debug(std::cout << "derivative mode\n");
        Bezier dbz = derivative(bz);
    
        debug(std::cout << "initial = " << dbz << std::endl);
        std::vector<double> dsolutions = dbz.roots(Interval(left_t, right_t));
        debug(std::cout << "dsolutions = " << dsolutions << std::endl);
        
        double dsplit_t = 0.5;
        if(!dsolutions.empty()) {
            dsplit_t = dsolutions[0];
            split_t = left_t + (right_t - left_t)*dsplit_t;
            debug(std::cout << "split_value = " << bz(split_t) << std::endl);
            debug(std::cout << "spliting around " << dsplit_t << " = " 
                  << split_t << "\n");
        
        }
        std::pair<Bezier, Bezier> LR = bz.subdivide(dsplit_t);
        Left = LR.first;
        Right = LR.second;
    } else {
        // split at midpoint, because it is cheap
        Left[0] = Right[0];
        for (size_t i = 1; i < bz.size(); ++i)
        {
            for (size_t j = 0; j < bz.size()-i; ++j)
            {
                Right[j] = (Right[j] + Right[j+1]) * 0.5;
            }
            Left[i] = Right[0];
        }
    }
    debug(std::cout << "Solution is exactly on the subdivision point.\n");
    debug(std::cout << Left << " , " << Right << std::endl);
    Left = reverse(Left);
    while(Right.order() > 0 and fabs(Right[0]) <= 1e-10) {
        debug(std::cout << "deflate\n");
        Right = Right.deflate();
        Left = Left.deflate();
        solutions.push_back(split_t);
    }
    Left = reverse(Left);
    if (Right.order() > 0) {
        debug(std::cout << Left << " , " << Right << std::endl);
        find_bernstein_roots(Left, depth+1, left_t, split_t);
        find_bernstein_roots(Right, depth+1, split_t, right_t);
    }
}