1 /* From Sanchez-Reyes 1997
2 W_{j,k} = W_{n0j, n-k} = choose(n-2k-1, j-k)choose(2k+1,k)/choose(n,j)
3 for k=0,...,q-1; j = k, ...,n-k-1
4 W_{q,q} = 1 (n even)
6 This is wrong, it should read
7 W_{j,k} = W_{n0j, n-k} = choose(n-2k-1, j-k)/choose(n,j)
8 for k=0,...,q-1; j = k, ...,n-k-1
9 W_{q,q} = 1 (n even)
11 */
12 #include <2geom/sbasis-to-bezier.h>
13 #include <2geom/choose.h>
14 #include <2geom/svg-path.h>
15 #include <iostream>
16 #include <2geom/exception.h>
18 namespace Geom{
20 double W(unsigned n, unsigned j, unsigned k) {
21 unsigned q = (n+1)/2;
22 if((n & 1) == 0 && j == q && k == q)
23 return 1;
24 if(k > n-k) return W(n, n-j, n-k);
25 assert((k <= q));
26 if(k >= q) return 0;
27 //assert(!(j >= n-k));
28 if(j >= n-k) return 0;
29 //assert(!(j < k));
30 if(j < k) return 0;
31 return choose<double>(n-2*k-1, j-k) /
32 choose<double>(n,j);
33 }
35 // this produces a degree 2q bezier from a degree k sbasis
36 Bezier
37 sbasis_to_bezier(SBasis const &B, unsigned q) {
38 if(q == 0) {
39 q = B.size();
40 /*if(B.back()[0] == B.back()[1]) {
41 n--;
42 }*/
43 }
44 unsigned n = q*2;
45 Bezier result = Bezier(Bezier::Order(n-1));
46 if(q > B.size())
47 q = B.size();
48 n--;
49 for(unsigned k = 0; k < q; k++) {
50 for(unsigned j = 0; j <= n-k; j++) {
51 result[j] += (W(n, j, k)*B[k][0] +
52 W(n, n-j, k)*B[k][1]);
53 }
54 }
55 return result;
56 }
58 double mopi(int i) {
59 return (i&1)?-1:1;
60 }
62 // this produces a degree k sbasis from a degree 2q bezier
63 SBasis
64 bezier_to_sbasis(Bezier const &B) {
65 unsigned n = B.size();
66 unsigned q = (n+1)/2;
67 SBasis result;
68 result.resize(q+1);
69 for(unsigned k = 0; k < q; k++) {
70 result[k][0] = result[k][1] = 0;
71 for(unsigned j = 0; j <= n-k; j++) {
72 result[k][0] += mopi(int(j)-int(k))*W(n, j, k)*B[j];
73 result[k][1] += mopi(int(j)-int(k))*W(n, j, k)*B[j];
74 //W(n, n-j, k)*B[k][1]);
75 }
76 }
77 return result;
78 }
80 // this produces a 2q point bezier from a degree q sbasis
81 std::vector<Geom::Point>
82 sbasis_to_bezier(D2<SBasis> const &B, unsigned qq) {
83 std::vector<Geom::Point> result;
84 if(qq == 0) {
85 qq = sbasis_size(B);
86 }
87 unsigned n = qq * 2;
88 result.resize(n, Geom::Point(0,0));
89 n--;
90 for(unsigned dim = 0; dim < 2; dim++) {
91 unsigned q = qq;
92 if(q > B[dim].size())
93 q = B[dim].size();
94 for(unsigned k = 0; k < q; k++) {
95 for(unsigned j = 0; j <= n-k; j++) {
96 result[j][dim] += (W(n, j, k)*B[dim][k][0] +
97 W(n, n-j, k)*B[dim][k][1]);
98 }
99 }
100 }
101 return result;
102 }
103 /*
104 template <unsigned order>
105 D2<Bezier<order> > sbasis_to_bezier(D2<SBasis> const &B) {
106 return D2<Bezier<order> >(sbasis_to_bezier<order>(B[0]), sbasis_to_bezier<order>(B[1]));
107 }
108 */
110 #if 0 // using old path
111 //std::vector<Geom::Point>
112 // mutating
113 void
114 subpath_from_sbasis(Geom::OldPathSetBuilder &pb, D2<SBasis> const &B, double tol, bool initial) {
115 assert(B.IS_FINITE());
116 if(B.tail_error(2) < tol || B.size() == 2) { // nearly cubic enough
117 if(B.size() == 1) {
118 if (initial) {
119 pb.start_subpath(Geom::Point(B[0][0][0], B[1][0][0]));
120 }
121 pb.push_line(Geom::Point(B[0][0][1], B[1][0][1]));
122 } else {
123 std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
124 if (initial) {
125 pb.start_subpath(bez[0]);
126 }
127 pb.push_cubic(bez[1], bez[2], bez[3]);
128 }
129 } else {
130 subpath_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol, initial);
131 subpath_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol, false);
132 }
133 }
135 /*
136 * This version works by inverting a reasonable upper bound on the error term after subdividing the
137 * curve at $a$. We keep biting off pieces until there is no more curve left.
138 *
139 * Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$. A
140 * subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$. Let this be the desired
141 * tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$
142 */
143 void
144 subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, double tol, bool initial) {
145 const unsigned k = 2; // cubic bezier
146 double te = B.tail_error(k);
147 assert(B[0].IS_FINITE());
148 assert(B[1].IS_FINITE());
150 //std::cout << "tol = " << tol << std::endl;
151 while(1) {
152 double A = std::sqrt(tol/te); // pow(te, 1./k)
153 double a = A;
154 if(A < 1) {
155 A = std::min(A, 0.25);
156 a = 0.5 - std::sqrt(0.25 - A); // quadratic formula
157 if(a > 1) a = 1; // clamp to the end of the segment
158 } else
159 a = 1;
160 assert(a > 0);
161 //std::cout << "te = " << te << std::endl;
162 //std::cout << "A = " << A << "; a=" << a << std::endl;
163 D2<SBasis> Bs = compose(B, Linear(0, a));
164 assert(Bs.tail_error(k));
165 std::vector<Geom::Point> bez = sbasis_to_bezier(Bs, 2);
166 reverse(bez.begin(), bez.end());
167 if (initial) {
168 pb.start_subpath(bez[0]);
169 initial = false;
170 }
171 pb.push_cubic(bez[1], bez[2], bez[3]);
173 // move to next piece of curve
174 if(a >= 1) break;
175 B = compose(B, Linear(a, 1));
176 te = B.tail_error(k);
177 }
178 }
180 #endif
182 /*
183 * If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
184 */
185 void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
186 if (!B.isFinite()) {
187 THROW_EXCEPTION("assertion failed: B.isFinite()");
188 }
189 if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
190 if( !only_cubicbeziers && (sbasis_size(B) <= 1) ) {
191 pb.lineTo(B.at1());
192 } else {
193 std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
194 pb.curveTo(bez[1], bez[2], bez[3]);
195 }
196 } else {
197 build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol, only_cubicbeziers);
198 build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol, only_cubicbeziers);
199 }
200 }
202 /*
203 * If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
204 */
205 Path
206 path_from_sbasis(D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
207 PathBuilder pb;
208 pb.moveTo(B.at0());
209 build_from_sbasis(pb, B, tol, only_cubicbeziers);
210 pb.finish();
211 return pb.peek().front();
212 }
214 /*
215 * If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
216 */
217 //TODO: some of this logic should be lifted into svg-path
218 std::vector<Geom::Path>
219 path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol, bool only_cubicbeziers) {
220 Geom::PathBuilder pb;
221 if(B.size() == 0) return pb.peek();
222 Geom::Point start = B[0].at0();
223 pb.moveTo(start);
224 for(unsigned i = 0; ; i++) {
225 if(i+1 == B.size() || !are_near(B[i+1].at0(), B[i].at1(), tol)) {
226 //start of a new path
227 if(are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
228 pb.closePath();
229 //last line seg already there (because of .closePath())
230 goto no_add;
231 }
232 build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
233 if(are_near(start, B[i].at1())) {
234 //it's closed, the last closing segment was not a straight line so it needed to be added, but still make it closed here with degenerate straight line.
235 pb.closePath();
236 }
237 no_add:
238 if(i+1 >= B.size()) break;
239 start = B[i+1].at0();
240 pb.moveTo(start);
241 } else {
242 build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
243 }
244 }
245 pb.finish();
246 return pb.peek();
247 }
249 };
251 /*
252 Local Variables:
253 mode:c++
254 c-file-style:"stroustrup"
255 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
256 indent-tabs-mode:nil
257 fill-column:99
258 End:
259 */
260 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :