GCC Code Coverage Report

Directory: ./
File: src/2geom/solve-bezier.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 91 91 100.0%
Functions: 4 4 100.0%
Branches: 93 136 68.4%
Line Branch Exec Source
1 #include <2geom/solver.h>
2 #include <2geom/choose.h>
3 #include <2geom/bezier.h>
4 #include <2geom/point.h>
5
6 #include <cmath>
7 #include <algorithm>
8
9 /*** Find the zeros of a Bezier. The code subdivides until it is happy with the linearity of the
10 * function. This requires an O(degree^2) subdivision for each step, even when there is only one
11 * solution.
12 *
13 * We try fairly hard to correctly handle multiple roots.
14 */
15
16 //#define debug(x) do{x;}while(0)
17 #define debug(x)
18
19 namespace Geom {
20 namespace {
21
22 struct Bernsteins
23 {
24 static constexpr size_t MAX_DEPTH = 22;
25 std::vector<double> &solutions;
26
27 65241 Bernsteins(std::vector<double> &sol)
28 65241 : solutions{sol}
29 65241 {}
30
31 double secant(Bezier const &bz);
32
33 void find_bernstein_roots(Bezier const &bz, unsigned depth,
34 double left_t, double right_t);
35 };
36
37 } // namespace
38
39 void
40 67906 Bezier::find_bezier_roots(std::vector<double> &solutions,
41 double left_t, double right_t) const {
42
1/2
✓ Branch 2 taken 67906 times.
✗ Branch 3 not taken.
 67906 Bezier bz = *this;
43
44 // a constant bezier, even if identically zero, has no roots
45
3/4
✓ Branch 1 taken 67906 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 241 times.
✓ Branch 4 taken 67665 times.
 67906 if (bz.isConstant()) {
46 241 return;
47 }
48
49
2/2
✓ Branch 1 taken 5170 times.
✓ Branch 2 taken 67665 times.
 72835 while(bz[0] == 0) {
50 debug(std::cout << "deflate\n");
51
2/4
✓ Branch 2 taken 5170 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5170 times.
✗ Branch 6 not taken.
 5170 bz = bz.deflate();
52
1/2
✓ Branch 2 taken 5170 times.
✗ Branch 3 not taken.
 5170 solutions.push_back(0);
53 }
54
3/4
✓ Branch 1 taken 67665 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2424 times.
✓ Branch 4 taken 65241 times.
 67665 if (bz.degree() == 1) {
55 debug(std::cout << "linear\n");
56
57
2/2
✓ Branch 4 taken 2348 times.
✓ Branch 5 taken 76 times.
 2424 if (Geom::sgn(bz[0]) != Geom::sgn(bz[1])) {
58 2348 double d = bz[0] - bz[1];
59
1/2
✓ Branch 0 taken 2348 times.
✗ Branch 1 not taken.
 2348 if(d != 0) {
60 2348 double r = bz[0] / d;
61
2/4
✓ Branch 0 taken 2348 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2348 times.
✗ Branch 3 not taken.
 2348 if(0 <= r && r <= 1)
62
1/2
✓ Branch 1 taken 2348 times.
✗ Branch 2 not taken.
 2348 solutions.push_back(r);
63 }
64 }
65 2424 return;
66 }
67
68 //std::cout << "initial = " << bz << std::endl;
69 65241 Bernsteins B(solutions);
70
1/2
✓ Branch 1 taken 65241 times.
✗ Branch 2 not taken.
 65241 B.find_bernstein_roots(bz, 0, left_t, right_t);
71 //std::cout << solutions << std::endl;
72
2/2
✓ Branch 1 taken 65241 times.
✓ Branch 2 taken 2665 times.
 67906 }
73
74 172835 void Bernsteins::find_bernstein_roots(Bezier const &bz,
75 unsigned depth,
76 double left_t,
77 double right_t)
78 {
79 debug(std::cout << left_t << ", " << right_t << std::endl);
80 172835 size_t n_crossings = 0;
81
82 172835 int old_sign = Geom::sgn(bz[0]);
83 //std::cout << "w[0] = " << bz[0] << std::endl;
84
3/4
✓ Branch 1 taken 1013508 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 840673 times.
✓ Branch 4 taken 172835 times.
 1013508 for (size_t i = 1; i < bz.size(); i++)
85 {
86 //std::cout << "w[" << i << "] = " << w[i] << std::endl;
87 840673 int sign = Geom::sgn(bz[i]);
88
2/2
✓ Branch 0 taken 833317 times.
✓ Branch 1 taken 7356 times.
 840673 if (sign != 0)
89 {
90
3/4
✓ Branch 0 taken 233263 times.
✓ Branch 1 taken 600054 times.
✓ Branch 2 taken 233263 times.
✗ Branch 3 not taken.
 833317 if (sign != old_sign && old_sign != 0)
91 {
92 233263 ++n_crossings;
93 }
94 833317 old_sign = sign;
95 }
96 }
97 // if last control point is zero, that counts as crossing too
98
3/4
✓ Branch 1 taken 172835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7349 times.
✓ Branch 5 taken 165486 times.
 172835 if (bz[bz.size() - 1] == 0) {
99 7349 ++n_crossings;
100 }
101
102 //std::cout << "n_crossings = " << n_crossings << std::endl;
103
2/2
✓ Branch 0 taken 66205 times.
✓ Branch 1 taken 106630 times.
 172835 if (n_crossings == 0) return; // no solutions here
104
105
2/2
✓ Branch 0 taken 52655 times.
✓ Branch 1 taken 53975 times.
 106630 if (n_crossings == 1) /* Unique solution */
106 {
107 //std::cout << "depth = " << depth << std::endl;
108 /* Stop recursion when the tree is deep enough */
109 /* if deep enough, return 1 solution at midpoint */
110
2/2
✓ Branch 0 taken 72 times.
✓ Branch 1 taken 52583 times.
 52655 if (depth > MAX_DEPTH)
111 {
112 //printf("bottom out %d\n", depth);
113 72 const double Ax = right_t - left_t;
114
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 const double Ay = bz.at1() - bz.at0();
115
116
1/2
✓ Branch 3 taken 72 times.
✗ Branch 4 not taken.
 72 solutions.push_back(left_t - Ax*bz.at0() / Ay);
117 72 return;
118 }
119
120
1/2
✓ Branch 1 taken 52583 times.
✗ Branch 2 not taken.
 52583 double r = secant(bz);
121
1/2
✓ Branch 2 taken 52583 times.
✗ Branch 3 not taken.
 52583 solutions.push_back(r*right_t + (1-r)*left_t);
122 52583 return;
123 }
124 /* Otherwise, solve recursively after subdividing control polygon */
125
1/2
✓ Branch 2 taken 53975 times.
✗ Branch 3 not taken.
 53975 Bezier::Order o(bz);
126
2/4
✓ Branch 2 taken 53975 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 53975 times.
✗ Branch 7 not taken.
 53975 Bezier Left(o), Right = bz;
127 53975 double split_t = (left_t + right_t) * 0.5;
128
129 // If subdivision is working poorly, split around the leftmost root of the derivative
130
2/2
✓ Branch 0 taken 22188 times.
✓ Branch 1 taken 31787 times.
 53975 if (depth > 2) {
131 debug(std::cout << "derivative mode\n");
132
1/2
✓ Branch 2 taken 22188 times.
✗ Branch 3 not taken.
 22188 Bezier dbz = derivative(bz);
133
134 debug(std::cout << "initial = " << dbz << std::endl);
135
2/4
✓ Branch 3 taken 22188 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22188 times.
✗ Branch 7 not taken.
 22188 std::vector<double> dsolutions = dbz.roots(Interval(left_t, right_t));
136 debug(std::cout << "dsolutions = " << dsolutions << std::endl);
137
138 22188 double dsplit_t = 0.5;
139
2/2
✓ Branch 1 taken 18291 times.
✓ Branch 2 taken 3897 times.
 22188 if(!dsolutions.empty()) {
140 18291 dsplit_t = dsolutions[0];
141 18291 split_t = left_t + (right_t - left_t)*dsplit_t;
142 debug(std::cout << "split_value = " << bz(split_t) << std::endl);
143 debug(std::cout << "splitting around " << dsplit_t << " = "
144 << split_t << "\n");
145
146 }
147
1/2
✓ Branch 2 taken 22188 times.
✗ Branch 3 not taken.
 22188 std::pair<Bezier, Bezier> LR = bz.subdivide(dsplit_t);
148
1/2
✓ Branch 1 taken 22188 times.
✗ Branch 2 not taken.
 22188 Left = LR.first;
149
1/2
✓ Branch 1 taken 22188 times.
✗ Branch 2 not taken.
 22188 Right = LR.second;
150 22188 } else {
151 // split at midpoint, because it is cheap
152 31787 Left[0] = Right[0];
153
3/4
✓ Branch 1 taken 193243 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 161456 times.
✓ Branch 4 taken 31787 times.
 193243 for (size_t i = 1; i < bz.size(); ++i)
154 {
155
3/4
✓ Branch 1 taken 911933 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 750477 times.
✓ Branch 4 taken 161456 times.
 911933 for (size_t j = 0; j < bz.size()-i; ++j)
156 {
157 750477 Right[j] = (Right[j] + Right[j+1]) * 0.5;
158 }
159 161456 Left[i] = Right[0];
160 }
161 }
162 debug(std::cout << "Solution is exactly on the subdivision point.\n");
163 debug(std::cout << Left << " , " << Right << std::endl);
164
2/4
✓ Branch 2 taken 53975 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 53975 times.
✗ Branch 6 not taken.
 53975 Left = reverse(Left);
165
7/8
✓ Branch 1 taken 71556 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 71378 times.
✓ Branch 4 taken 178 times.
✓ Branch 6 taken 17581 times.
✓ Branch 7 taken 53797 times.
✓ Branch 8 taken 17581 times.
✓ Branch 9 taken 53975 times.
 71556 while(Right.order() > 0 && fabs(Right[0]) <= 1e-10) {
166 debug(std::cout << "deflate\n");
167
2/4
✓ Branch 2 taken 17581 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 17581 times.
✗ Branch 6 not taken.
 17581 Right = Right.deflate();
168
2/4
✓ Branch 2 taken 17581 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 17581 times.
✗ Branch 6 not taken.
 17581 Left = Left.deflate();
169
1/2
✓ Branch 1 taken 17581 times.
✗ Branch 2 not taken.
 17581 solutions.push_back(split_t);
170 }
171
2/4
✓ Branch 2 taken 53975 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 53975 times.
✗ Branch 6 not taken.
 53975 Left = reverse(Left);
172
3/4
✓ Branch 1 taken 53975 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 53797 times.
✓ Branch 4 taken 178 times.
 53975 if (Right.order() > 0) {
173 debug(std::cout << Left << " , " << Right << std::endl);
174
1/2
✓ Branch 1 taken 53797 times.
✗ Branch 2 not taken.
 53797 find_bernstein_roots(Left, depth+1, left_t, split_t);
175
1/2
✓ Branch 1 taken 53797 times.
✗ Branch 2 not taken.
 53797 find_bernstein_roots(Right, depth+1, split_t, right_t);
176 }
177 53975 }
178
179 52583 double Bernsteins::secant(Bezier const &bz) {
180 52583 double s = 0, t = 1;
181 52583 double e = 1e-14;
182 52583 int side = 0;
183 52583 double r, fs = bz.at0(), ft = bz.at1();
184
185
1/2
✓ Branch 0 taken 577795 times.
✗ Branch 1 not taken.
 577795 for (size_t n = 0; n < 100; ++n)
186 {
187 577795 r = (fs*t - ft*s) / (fs - ft);
188
2/2
✓ Branch 0 taken 41986 times.
✓ Branch 1 taken 535809 times.
 577795 if (fabs(t-s) < e * fabs(t+s)) {
189 debug(std::cout << "error small " << fabs(t-s)
190 << ", accepting solution " << r
191 << "after " << n << "iterations\n");
192 41986 return r;
193 }
194
195 535809 double fr = bz.valueAt(r);
196
197
2/2
✓ Branch 0 taken 316210 times.
✓ Branch 1 taken 219599 times.
 535809 if (fr * ft > 0)
198 {
199 316210 t = r; ft = fr;
200
2/2
✓ Branch 0 taken 148322 times.
✓ Branch 1 taken 167888 times.
 316210 if (side == -1) fs /= 2;
201 316210 side = -1;
202 }
203
2/2
✓ Branch 0 taken 209002 times.
✓ Branch 1 taken 10597 times.
 219599 else if (fs * fr > 0)
204 {
205 209002 s = r; fs = fr;
206
2/2
✓ Branch 0 taken 45498 times.
✓ Branch 1 taken 163504 times.
 209002 if (side == +1) ft /= 2;
207 209002 side = +1;
208 }
209 10597 else break;
210 }
211 10597 return r;
212 }
213
214 } // namespace Geom
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:fileencoding=utf-8:textwidth=99 :
226