src/2geom/solve-bezier-parametric.cpp

   1 #include <2geom/solver.h>
   2 #include <2geom/point.h>
   3 #include <algorithm>
   4
   5 namespace  Geom{
   6
   7 /*** Find the zeros of the parametric function in 2d defined by two beziers X(t), Y(t).  The code subdivides until it happy with the linearity of the bezier.  This requires an n^2 subdivision for each step, even when there is only one solution.
   8  *
   9  * Perhaps it would be better to subdivide particularly around nodes with changing sign, rather than simply cutting in half.
  10  */
  11
  12 #define SGN(a)      (((a)<0) ? -1 : 1)
  13
  14 /*
  15  *  Forward declarations
  16  */
  17 static Geom::Point
  18 Bezier(Geom::Point const *V,
  19        unsigned degree,
  20        double t,
  21        Geom::Point *Left,
  22        Geom::Point *Right);
  23
  24 unsigned
  25 crossing_count(Geom::Point const *V, unsigned degree);
  26 static unsigned
  27 control_poly_flat_enough(Geom::Point const *V, unsigned degree);
  28 static double
  29 compute_x_intercept(Geom::Point const *V, unsigned degree);
  30
  31 const unsigned MAXDEPTH = 64;   /*  Maximum depth for recursion */
  32
  33 const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
  34
  35 unsigned total_steps, total_subs;
  36
  37 /*
  38  *  find_bezier_roots : Given an equation in Bernstein-Bezier form, find all
  39  *    of the roots in the interval [0, 1].  Return the number of roots found.
  40  */
  41 void
  42 find_parametric_bezier_roots(Geom::Point const *w, /* The control points  */
  43                   unsigned degree,      /* The degree of the polynomial */
  44                   std::vector<double> &solutions, /* RETURN candidate t-values */
  45                   unsigned depth)       /* The depth of the recursion */
  46 {
  47     total_steps++;
  48     const unsigned max_crossings = crossing_count(w, degree);
  49     switch (max_crossings) {
  50     case 0:     /* No solutions here    */
  51         return;
  52
  53     case 1:
  54         /* Unique solution      */
  55         /* Stop recursion when the tree is deep enough  */
  56         /* if deep enough, return 1 solution at midpoint  */
  57         if (depth >= MAXDEPTH) {
  58             solutions.push_back((w[0][Geom::X] + w[degree][Geom::X]) / 2.0);
  59             return;
  60         }
  61
  62         // I thought secant method would be faster here, but it'aint. -- njh
  63
  64         if (control_poly_flat_enough(w, degree)) {
  65             solutions.push_back(compute_x_intercept(w, degree));
  66             return;
  67         }
  68         break;
  69     }
  70
  71     /* Otherwise, solve recursively after subdividing control polygon  */
  72     Geom::Point Left[degree+1], /* New left and right  */
  73         Right[degree+1];        /* control polygons  */
  74     Bezier(w, degree, 0.5, Left, Right);
  75     total_subs ++;
  76     find_parametric_bezier_roots(Left,  degree, solutions, depth+1);
  77     find_parametric_bezier_roots(Right, degree, solutions, depth+1);
  78 }
  79
  80
  81 /*
  82  * crossing_count:
  83  *  Count the number of times a Bezier control polygon
  84  *  crosses the 0-axis. This number is >= the number of roots.
  85  *
  86  */
  87 unsigned
  88 crossing_count(Geom::Point const *V,    /*  Control pts of Bezier curve */
  89                unsigned degree) /*  Degree of Bezier curve      */
  90 {
  91     unsigned    n_crossings = 0;        /*  Number of zero-crossings */
  92
  93     int old_sign = SGN(V[0][Geom::Y]);
  94     for (unsigned i = 1; i <= degree; i++) {
  95         int sign = SGN(V[i][Geom::Y]);
  96         if (sign != old_sign)
  97             n_crossings++;
  98         old_sign = sign;
  99     }
 100     return n_crossings;
 101 }
 102
 103
 104
 105 /*
 106  *  control_poly_flat_enough :
 107  *      Check if the control polygon of a Bezier curve is flat enough
 108  *      for recursive subdivision to bottom out.
 109  *
 110  */
 111 static unsigned
 112 control_poly_flat_enough(Geom::Point const *V, /* Control points        */
 113                          unsigned degree)       /* Degree of polynomial */
 114 {
 115     /* Find the perpendicular distance from each interior control point to line connecting V[0] and
 116      * V[degree] */
 117
 118     /* Derive the implicit equation for line connecting first */
 119     /*  and last control points */
 120     const double a = V[0][Geom::Y] - V[degree][Geom::Y];
 121     const double b = V[degree][Geom::X] - V[0][Geom::X];
 122     const double c = V[0][Geom::X] * V[degree][Geom::Y] - V[degree][Geom::X] * V[0][Geom::Y];
 123
 124     const double abSquared = (a * a) + (b * b);
 125
 126     double distance[degree]; /* Distances from pts to line */
 127     for (unsigned i = 1; i < degree; i++) {
 128         /* Compute distance from each of the points to that line */
 129         double & dist(distance[i-1]);
 130         const double d = a * V[i][Geom::X] + b * V[i][Geom::Y] + c;
 131         dist = d*d / abSquared;
 132         if (d < 0.0)
 133             dist = -dist;
 134     }
 135
 136
 137     // Find the largest distance
 138     double max_distance_above = 0.0;
 139     double max_distance_below = 0.0;
 140     for (unsigned i = 0; i < degree-1; i++) {
 141         const double d = distance[i];
 142         if (d < 0.0)
 143             max_distance_below = std::min(max_distance_below, d);
 144         if (d > 0.0)
 145             max_distance_above = std::max(max_distance_above, d);
 146     }
 147
 148     const double intercept_1 = (c + max_distance_above) / -a;
 149     const double intercept_2 = (c + max_distance_below) / -a;
 150
 151     /* Compute bounding interval*/
 152     const double left_intercept = std::min(intercept_1, intercept_2);
 153     const double right_intercept = std::max(intercept_1, intercept_2);
 154
 155     const double error = 0.5 * (right_intercept - left_intercept);
 156
 157     if (error < BEPSILON)
 158         return 1;
 159
 160     return 0;
 161 }
 162
 163
 164
 165 /*
 166  *  compute_x_intercept :
 167  *      Compute intersection of chord from first control point to last
 168  *      with 0-axis.
 169  *
 170  */
 171 static double
 172 compute_x_intercept(Geom::Point const *V, /*  Control points    */
 173                     unsigned degree) /*  Degree of curve        */
 174 {
 175     const Geom::Point A = V[degree] - V[0];
 176
 177     return (A[Geom::X]*V[0][Geom::Y] - A[Geom::Y]*V[0][Geom::X]) / -A[Geom::Y];
 178 }
 179
 180
 181 /*
 182  *  Bezier :
 183  *      Evaluate a Bezier curve at a particular parameter value
 184  *      Fill in control points for resulting sub-curves.
 185  *
 186  */
 187 static Geom::Point
 188 Bezier(Geom::Point const *V, /* Control pts     */
 189        unsigned degree, /* Degree of bezier curve */
 190        double t,        /* Parameter value */
 191        Geom::Point *Left,       /* RETURN left half ctl pts */
 192        Geom::Point *Right)      /* RETURN right half ctl pts */
 193 {
 194     Geom::Point Vtemp[degree+1][degree+1];
 195
 196     /* Copy control points      */
 197     std::copy(V, V+degree+1, Vtemp[0]);
 198
 199     /* Triangle computation     */
 200     for (unsigned i = 1; i <= degree; i++) {
 201         for (unsigned j = 0; j <= degree - i; j++) {
 202             Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
 203         }
 204     }
 205
 206     for (unsigned j = 0; j <= degree; j++)
 207         Left[j]  = Vtemp[j][0];
 208     for (unsigned j = 0; j <= degree; j++)
 209         Right[j] = Vtemp[degree-j][j];
 210
 211     return (Vtemp[degree][0]);
 212 }
 213
 214 };
 215
 216 /*
 217   Local Variables:
 218   mode:c++
 219   c-file-style:"stroustrup"
 220   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 221   indent-tabs-mode:nil
 222   fill-column:99
 223   End:
 224 */
 225 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :