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

   1 #include "solver.h"
   2 #include "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 /*
  16  *  Forward declarations
  17  */
  18 static void
  19 Bernstein(double const *V,
  20           unsigned degree,
  21           double t,
  22           double *Left,
  23           double *Right);
  24
  25 static unsigned
  26 control_poly_flat_enough(double const *V, unsigned degree,
  27                          double left_t, double right_t);
  28
  29 const unsigned MAXDEPTH = 64;   /*  Maximum depth for recursion */
  30
  31 const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
  32
  33 /*
  34  *  find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all
  35  *    of the roots in the open interval (0, 1).  Return the number of roots found.
  36  */
  37 void
  38 find_bernstein_roots(double const *w, /* The control points  */
  39                      unsigned degree,   /* The degree of the polynomial */
  40                      std::vector<double> &solutions, /* RETURN candidate t-values */
  41                      unsigned depth,    /* The depth of the recursion */
  42                      double left_t, double right_t)
  43 {
  44     unsigned    n_crossings = 0;        /*  Number of zero-crossings */
  45
  46     int old_sign = SGN(w[0]);
  47     for (unsigned i = 1; i <= degree; i++) {
  48         int sign = SGN(w[i]);
  49         if (sign) {
  50             if (sign != old_sign && old_sign) {
  51                n_crossings++;
  52             }
  53             old_sign = sign;
  54         }
  55     }
  56
  57     switch (n_crossings) {
  58     case 0:     /* No solutions here    */
  59         return;
  60
  61     case 1:
  62         /* Unique solution      */
  63         /* Stop recursion when the tree is deep enough  */
  64         /* if deep enough, return 1 solution at midpoint  */
  65         if (depth >= MAXDEPTH) {
  66             solutions.push_back((left_t + right_t) / 2.0);
  67             return;
  68         }
  69
  70         // I thought secant method would be faster here, but it'aint. -- njh
  71
  72         if (control_poly_flat_enough(w, degree, left_t, right_t)) {
  73             const double Ax = right_t - left_t;
  74             const double Ay = w[degree] - w[0];
  75
  76             solutions.push_back(left_t - Ax*w[0] / Ay);
  77             return;
  78         }
  79         break;
  80     }
  81
  82     /* Otherwise, solve recursively after subdividing control polygon  */
  83     double Left[degree+1],      /* New left and right  */
  84            Right[degree+1];     /* control polygons  */
  85     const double split = 0.5;
  86     Bernstein(w, degree, split, Left, Right);
  87
  88     double mid_t = left_t*(1-split) + right_t*split;
  89
  90     find_bernstein_roots(Left,  degree, solutions, depth+1, left_t, mid_t);
  91
  92     /* Solution is exactly on the subdivision point. */
  93     if (Right[0] == 0)
  94         solutions.push_back(mid_t);
  95
  96     find_bernstein_roots(Right, degree, solutions, depth+1, mid_t, right_t);
  97 }
  98
  99 /*
 100  *  control_poly_flat_enough :
 101  *      Check if the control polygon of a Bernstein curve is flat enough
 102  *      for recursive subdivision to bottom out.
 103  *
 104  */
 105 static unsigned
 106 control_poly_flat_enough(double const *V, /* Control points     */
 107                          unsigned degree,
 108                          double left_t, double right_t) /* Degree of polynomial */
 109 {
 110     /* Find the perpendicular distance from each interior control point to line connecting V[0] and
 111      * V[degree] */
 112
 113     /* Derive the implicit equation for line connecting first */
 114     /*  and last control points */
 115     const double a = V[0] - V[degree];
 116     const double b = right_t - left_t;
 117     const double c = left_t * V[degree] - right_t * V[0] + a * left_t;
 118
 119     double max_distance_above = 0.0;
 120     double max_distance_below = 0.0;
 121     double ii = 0, dii = 1./degree;
 122     for (unsigned i = 1; i < degree; i++) {
 123         ii += dii;
 124         /* Compute distance from each of the points to that line */
 125         const double d = (a + V[i]) * ii*b  + c;
 126         double dist = d*d;
 127     // Find the largest distance
 128         if (d < 0.0)
 129             max_distance_below = std::min(max_distance_below, -dist);
 130         else
 131             max_distance_above = std::max(max_distance_above, dist);
 132     }
 133
 134     const double abSquared = (a * a) + (b * b);
 135
 136     const double intercept_1 = -(c + max_distance_above / abSquared);
 137     const double intercept_2 = -(c + max_distance_below / abSquared);
 138
 139     /* Compute bounding interval*/
 140     const double left_intercept = std::min(intercept_1, intercept_2);
 141     const double right_intercept = std::max(intercept_1, intercept_2);
 142
 143     const double error = 0.5 * (right_intercept - left_intercept);
 144
 145     if (error < BEPSILON * a)
 146         return 1;
 147
 148     return 0;
 149 }
 150
 151
 152
 153 /*
 154  *  Bernstein :
 155  *      Evaluate a Bernstein function at a particular parameter value
 156  *      Fill in control points for resulting sub-curves.
 157  *
 158  */
 159 static void
 160 Bernstein(double const *V, /* Control pts       */
 161           unsigned degree,      /* Degree of bernstein curve */
 162           double t,     /* Parameter value */
 163           double *Left, /* RETURN left half ctl pts */
 164           double *Right)        /* RETURN right half ctl pts */
 165 {
 166     double Vtemp[degree+1][degree+1];
 167
 168     /* Copy control points      */
 169     std::copy(V, V+degree+1, Vtemp[0]);
 170
 171     /* Triangle computation     */
 172     const double omt = (1-t);
 173     Left[0] = Vtemp[0][0];
 174     Right[degree] = Vtemp[0][degree];
 175     for (unsigned i = 1; i <= degree; i++) {
 176         for (unsigned j = 0; j <= degree - i; j++) {
 177             Vtemp[i][j] = omt*Vtemp[i-1][j] + t*Vtemp[i-1][j+1];
 178         }
 179         Left[i] = Vtemp[i][0];
 180         Right[degree-i] = Vtemp[i][degree-i];
 181     }
 182 }
 183
 184 };
 185
 186 /*
 187   Local Variables:
 188   mode:c++
 189   c-file-style:"stroustrup"
 190   c-file-offsets:((innamespace . 0)(substatement-open . 0))
 191   indent-tabs-mode:nil
 192   c-brace-offset:0
 193   fill-column:99
 194   End:
 195   vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
 196 */
 197