GCC Code Coverage Report

Directory: ./
File: src/2geom/solve-bezier-one-d.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 54 65 83.1%
Functions: 3 4 75.0%
Branches: 35 56 62.5%
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 the bernstein function. The code subdivides until it is happy with the
10 * linearity of the function. This requires an O(degree^2) subdivision for each step, even when
11 * there is only one solution.
12 */
13
14 namespace Geom {
15 namespace {
16
17 /**
18 * 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.
19 **/
20 struct Bernsteins
21 {
22 static constexpr size_t MAX_DEPTH = 53;
23 size_t degree, N;
24 std::vector<double> &solutions;
25
26 22188 Bernsteins(size_t _degree, std::vector<double> &sol)
27 22188 : degree{_degree}
28 22188 , N{degree + 1}
29 22188 , solutions{sol}
30 22188 {}
31
32 void
33 find_bernstein_roots(double const *w, /* The control points */
34 unsigned depth, /* The depth of the recursion */
35 double left_t, double right_t);
36 };
37
38 } // namespace
39
40 /*
41 * find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all
42 * of the roots in the open interval (0, 1). Return the number of roots found.
43 */
44 void
45 22188 find_bernstein_roots(double const *w, /* The control points */
46 unsigned degree, /* The degree of the polynomial */
47 std::vector<double> &solutions, /* RETURN candidate t-values */
48 unsigned depth, /* The depth of the recursion */
49 double left_t, double right_t, bool /*use_secant*/)
50 {
51 22188 Bernsteins B(degree, solutions);
52
1/2
✓ Branch 1 taken 22188 times.
✗ Branch 2 not taken.
 22188 B.find_bernstein_roots(w, depth, left_t, right_t);
53 22188 }
54
55 void
56 ✗ find_bernstein_roots(std::vector<double> &solutions, /* RETURN candidate t-values */
57 Geom::Bezier const &bz, /* The control points */
58 double left_t, double right_t)
59 {
60 ✗ Bernsteins B(bz.degree(), solutions);
61 ✗ Geom::Bezier& bzl = const_cast<Geom::Bezier&>(bz);
62 ✗ double* w = &(bzl[0]);
63 ✗ B.find_bernstein_roots(w, 0, left_t, right_t);
64 ✗ }
65
66 405298 void Bernsteins::find_bernstein_roots(double const *w, /* The control points */
67 unsigned depth, /* The depth of the recursion */
68 double left_t,
69 double right_t)
70 {
71 405298 size_t n_crossings = 0;
72
73 405298 int old_sign = Geom::sgn(w[0]);
74 //std::cout << "w[0] = " << w[0] << std::endl;
75
2/2
✓ Branch 0 taken 2819424 times.
✓ Branch 1 taken 405298 times.
 3224722 for (size_t i = 1; i < N; i++)
76 {
77 //std::cout << "w[" << i << "] = " << w[i] << std::endl;
78 2819424 int sign = Geom::sgn(w[i]);
79
1/2
✓ Branch 0 taken 2819424 times.
✗ Branch 1 not taken.
 2819424 if (sign != 0)
80 {
81
3/4
✓ Branch 0 taken 579462 times.
✓ Branch 1 taken 2239962 times.
✓ Branch 2 taken 579462 times.
✗ Branch 3 not taken.
 2819424 if (sign != old_sign && old_sign != 0)
82 {
83 579462 ++n_crossings;
84 }
85 2819424 old_sign = sign;
86 }
87 }
88 //std::cout << "n_crossings = " << n_crossings << std::endl;
89
2/2
✓ Branch 0 taken 193703 times.
✓ Branch 1 taken 211595 times.
 405298 if (n_crossings == 0) return; // no solutions here
90
91
2/2
✓ Branch 0 taken 20040 times.
✓ Branch 1 taken 191555 times.
 211595 if (n_crossings == 1) /* Unique solution */
92 {
93 //std::cout << "depth = " << depth << std::endl;
94 /* Stop recursion when the tree is deep enough */
95 /* if deep enough, return 1 solution at midpoint */
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 20040 times.
 20040 if (depth > MAX_DEPTH)
97 {
98 //printf("bottom out %d\n", depth);
99 ✗ const double Ax = right_t - left_t;
100 ✗ const double Ay = w[degree] - w[0];
101
102 ✗ solutions.push_back(left_t - Ax*w[0] / Ay);
103 ✗ return;
104 }
105
106
107 20040 double s = 0, t = 1;
108 20040 double e = 1e-10;
109 20040 int side = 0;
110 20040 double r, fs = w[0], ft = w[degree];
111
112
1/2
✓ Branch 0 taken 248344 times.
✗ Branch 1 not taken.
 248344 for (size_t n = 0; n < 100; ++n)
113 {
114 248344 r = (fs*t - ft*s) / (fs - ft);
115
2/2
✓ Branch 0 taken 19734 times.
✓ Branch 1 taken 228610 times.
 248344 if (fabs(t-s) < e * fabs(t+s)) break;
116
117 228610 double fr = bernstein_value_at(r, w, degree);
118
119
2/2
✓ Branch 0 taken 109154 times.
✓ Branch 1 taken 119456 times.
 228610 if (fr * ft > 0)
120 {
121 109154 t = r; ft = fr;
122
2/2
✓ Branch 0 taken 46063 times.
✓ Branch 1 taken 63091 times.
 109154 if (side == -1) fs /= 2;
123 109154 side = -1;
124 }
125
2/2
✓ Branch 0 taken 119150 times.
✓ Branch 1 taken 306 times.
 119456 else if (fs * fr > 0)
126 {
127 119150 s = r; fs = fr;
128
2/2
✓ Branch 0 taken 56877 times.
✓ Branch 1 taken 62273 times.
 119150 if (side == +1) ft /= 2;
129 119150 side = +1;
130 }
131 306 else break;
132 }
133
1/2
✓ Branch 2 taken 20040 times.
✗ Branch 3 not taken.
 20040 solutions.push_back(r*right_t + (1-r)*left_t);
134 20040 return;
135
136 }
137
138 /* Otherwise, solve recursively after subdividing control polygon */
139
2/4
✓ Branch 0 taken 191555 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 191555 times.
✗ Branch 5 not taken.
 191555 double* LR = new double[2*N];
140 191555 double* Left = LR;
141 191555 double* Right = LR + N;
142
143
1/2
✓ Branch 1 taken 191555 times.
✗ Branch 2 not taken.
 191555 std::copy(w, w + N, Right);
144
145 191555 Left[0] = Right[0];
146
2/2
✓ Branch 0 taken 1333376 times.
✓ Branch 1 taken 191555 times.
 1524931 for (size_t i = 1; i < N; ++i)
147 {
148
2/2
✓ Branch 0 taken 5402960 times.
✓ Branch 1 taken 1333376 times.
 6736336 for (size_t j = 0; j < N-i; ++j)
149 {
150 5402960 Right[j] = (Right[j] + Right[j+1]) * 0.5;
151 }
152 1333376 Left[i] = Right[0];
153 }
154
155 191555 double mid_t = (left_t + right_t) * 0.5;
156
157
158
1/2
✓ Branch 1 taken 191555 times.
✗ Branch 2 not taken.
 191555 find_bernstein_roots(Left, depth+1, left_t, mid_t);
159
160
161 /* Solution is exactly on the subdivision point. */
162
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 191555 times.
 191555 if (Right[0] == 0)
163 {
164 ✗ solutions.push_back(mid_t);
165 }
166
167
1/2
✓ Branch 1 taken 191555 times.
✗ Branch 2 not taken.
 191555 find_bernstein_roots(Right, depth+1, mid_t, right_t);
168
1/2
✓ Branch 0 taken 191555 times.
✗ Branch 1 not taken.
 191555 delete[] LR;
169 }
170
171 } // namespace Geom
172
173 /*
174 Local Variables:
175 mode:c++
176 c-file-style:"stroustrup"
177 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
178 indent-tabs-mode:nil
179 fill-column:99
180 End:
181 */
182 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
183