1 #include <2geom/solver.h>
2 #include <2geom/point.h>
3 #include <algorithm>
5 namespace Geom{
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 */
12 #define SGN(a) (((a)<0) ? -1 : 1)
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);
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);
31 const unsigned MAXDEPTH = 64; /* Maximum depth for recursion */
33 const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
35 unsigned total_steps, total_subs;
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;
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 }
62 // I thought secant method would be faster here, but it'aint. -- njh
64 if (control_poly_flat_enough(w, degree)) {
65 solutions.push_back(compute_x_intercept(w, degree));
66 return;
67 }
68 break;
69 }
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 }
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 */
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 }
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] */
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];
124 const double abSquared = (a * a) + (b * b);
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 }
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 }
148 const double intercept_1 = (c + max_distance_above) / -a;
149 const double intercept_2 = (c + max_distance_below) / -a;
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);
155 const double error = 0.5 * (right_intercept - left_intercept);
157 if (error < BEPSILON)
158 return 1;
160 return 0;
161 }
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];
177 return (A[Geom::X]*V[0][Geom::Y] - A[Geom::Y]*V[0][Geom::X]) / -A[Geom::Y];
178 }
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];
196 /* Copy control points */
197 std::copy(V, V+degree+1, Vtemp[0]);
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 }
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];
211 return (Vtemp[degree][0]);
212 }
214 };
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 :