src/2geom/solve-bezier-one-d.cpp

   1 #include <2geom/solver.h>
   2 #include <2geom/point.h>
   3 #include <algorithm>
   4
   5 /*** Find the zeros of the bernstein function.  The code subdivides until it is happy with the
   6  * linearity of the function.  This requires an O(degree^2) subdivision for each step, even when
   7  * there is only one solution.
   8  */
   9
  10 namespace Geom{
  11
  12 template<class t>
  13 static int SGN(t x) { return (x > 0 ? 1 : (x < 0 ? -1 : 0)); }
  14
  15 const unsigned MAXDEPTH = 23; // Maximum depth for recursion.  Using floats means 23 bits precision max
  16
  17 const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
  18
  19 /**
  20  * This function is called _a lot_.  We have included various manual memory management stuff to reduce the amount of mallocing that goes on.  In the future it is possible that this will hurt performance.
  21  **/
  22 class Bernsteins{
  23 public:
  24     double *Vtemp;
  25     unsigned N,degree;
  26     std::vector<double> &solutions;
  27     Bernsteins(int degr, std::vector<double> &so) : N(degr+1), degree(degr),solutions(so) {
  28         Vtemp = new double[N*2];
  29     }
  30     ~Bernsteins() {
  31         delete[] Vtemp;
  32     }
  33     void subdivide(double const *V,
  34                    double t,
  35                    double *Left,
  36                    double *Right);
  37
  38     unsigned
  39     control_poly_flat_enough(double const *V);
  40
  41     void
  42     find_bernstein_roots(double const *w, /* The control points  */
  43                          unsigned depth,  /* The depth of the recursion */
  44                          double left_t, double right_t);
  45 };
  46 /*
  47  *  find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all
  48  *    of the roots in the open interval (0, 1).  Return the number of roots found.
  49  */
  50 void
  51 find_bernstein_roots(double const *w, /* The control points  */
  52                      unsigned degree,   /* The degree of the polynomial */
  53                      std::vector<double> &solutions, /* RETURN candidate t-values */
  54                      unsigned depth,    /* The depth of the recursion */
  55                      double left_t, double right_t)
  56 {
  57     Bernsteins B(degree, solutions);
  58     B.find_bernstein_roots(w, depth, left_t, right_t);
  59 }
  60
  61 void
  62 Bernsteins::find_bernstein_roots(double const *w, /* The control points  */
  63                      unsigned depth,    /* The depth of the recursion */
  64                      double left_t, double right_t)
  65 {
  66
  67     unsigned    n_crossings = 0;        /*  Number of zero-crossings */
  68
  69     int old_sign = SGN(w[0]);
  70     for (unsigned i = 1; i < N; i++) {
  71         int sign = SGN(w[i]);
  72         if (sign) {
  73             if (sign != old_sign && old_sign) {
  74                n_crossings++;
  75             }
  76             old_sign = sign;
  77         }
  78     }
  79
  80     if (n_crossings == 0) // no solutions here
  81         return;
  82
  83     if (n_crossings == 1) {
  84         /* Unique solution      */
  85         /* Stop recursion when the tree is deep enough  */
  86         /* if deep enough, return 1 solution at midpoint  */
  87         if (depth >= MAXDEPTH) {
  88             //printf("bottom out %d\n", depth);
  89             const double Ax = right_t - left_t;
  90             const double Ay = w[degree] - w[0];
  91
  92             solutions.push_back(left_t - Ax*w[0] / Ay);
  93             return;
  94             solutions.push_back((left_t + right_t) / 2.0);
  95             return;
  96         }
  97
  98         // I thought secant method would be faster here, but it'aint. -- njh
  99
 100         // The original code went to some effort to try and terminate early when the curve is sufficiently flat.  However, it seems that in practice it almost always bottoms out, so leaving this code out actually speeds things up
 101         if(0) if (control_poly_flat_enough(w)) {
 102                 //printf("flatten out %d\n", depth);
 103                 const double Ax = right_t - left_t;
 104                 const double Ay = w[degree] - w[0];
 105
 106                 solutions.push_back(left_t - Ax*w[0] / Ay);
 107                 return;
 108             }
 109
 110     }
 111
 112     /* Otherwise, solve recursively after subdividing control polygon  */
 113     double Left[N],     /* New left and right  */
 114            Right[N];    /* control polygons  */
 115     const double t = 0.5;
 116
 117
 118 /*
 119  *  Bernstein :
 120  *      Evaluate a Bernstein function at a particular parameter value
 121  *      Fill in control points for resulting sub-curves.
 122  *
 123  */
 124     for (unsigned i = 0; i < N; i++)
 125         Vtemp[i] = w[i];
 126
 127     /* Triangle computation     */
 128     const double omt = (1-t);
 129     Left[0] = Vtemp[0];
 130     Right[degree] = Vtemp[degree];
 131     double *prev_row = Vtemp;
 132     double *row = Vtemp + N;
 133     for (unsigned i = 1; i < N; i++) {
 134         for (unsigned j = 0; j < N - i; j++) {
 135             row[j] = omt*prev_row[j] + t*prev_row[j+1];
 136         }
 137         Left[i] = row[0];
 138         Right[degree-i] = row[degree-i];
 139         std::swap(prev_row, row);
 140     }
 141
 142     double mid_t = left_t*(1-t) + right_t*t;
 143
 144     find_bernstein_roots(Left, depth+1, left_t, mid_t);
 145
 146     /* Solution is exactly on the subdivision point. */
 147     if (Right[0] == 0)
 148         solutions.push_back(mid_t);
 149
 150     find_bernstein_roots(Right, depth+1, mid_t, right_t);
 151 }
 152
 153 /*
 154  *  control_poly_flat_enough :
 155  *      Check if the control polygon of a Bernstein curve is flat enough
 156  *      for recursive subdivision to bottom out.
 157  *
 158  */
 159 unsigned
 160 Bernsteins::control_poly_flat_enough(double const *V)
 161 {
 162     /* Find the perpendicular distance from each interior control point to line connecting V[0] and
 163      * V[degree] */
 164
 165     /* Derive the implicit equation for line connecting first */
 166     /*  and last control points */
 167     const double a = V[0] - V[degree];
 168
 169     double max_distance_above = 0.0;
 170     double max_distance_below = 0.0;
 171     double ii = 0, dii = 1./degree;
 172     for (unsigned i = 1; i < degree; i++) {
 173         ii += dii;
 174         /* Compute distance from each of the points to that line */
 175         const double d = (a + V[i]) * ii  - a;
 176         double dist = d*d;
 177     // Find the largest distance
 178         if (d < 0.0)
 179             max_distance_below = std::min(max_distance_below, -dist);
 180         else
 181             max_distance_above = std::max(max_distance_above, dist);
 182     }
 183
 184     const double abSquared = 1./((a * a) + 1);
 185
 186     const double intercept_1 = (a - max_distance_above * abSquared);
 187     const double intercept_2 = (a - max_distance_below * abSquared);
 188
 189     /* Compute bounding interval*/
 190     const double left_intercept = std::min(intercept_1, intercept_2);
 191     const double right_intercept = std::max(intercept_1, intercept_2);
 192
 193     const double error = 0.5 * (right_intercept - left_intercept);
 194     //printf("error %g %g %g\n", error, a, BEPSILON * a);
 195     return error < BEPSILON * a;
 196 }
 197
 198
 199
 200 };
 201
 202 /*
 203   Local Variables:
 204   mode:c++
 205   c-file-style:"stroustrup"
 206   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 207   indent-tabs-mode:nil
 208   fill-column:99
 209   End:
 210 */
 211 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :