497091d1c65f29e7b95b4d4dc9d4e7ffd61ddd6d
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 "sbasis-to-bezier.h"
13 #include "choose.h"
14 #include "svg-path.h"
15 #include <iostream>
16 #include "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 void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol) {
183 if (!B.isFinite()) {
184 throwException("assertion failed: B.isFinite()");
185 }
186 if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
187 if(sbasis_size(B) <= 1) {
188 pb.lineTo(B.at1());
189 } else {
190 std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
191 pb.curveTo(bez[1], bez[2], bez[3]);
192 }
193 } else {
194 build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol);
195 build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol);
196 }
197 }
199 Path
200 path_from_sbasis(D2<SBasis> const &B, double tol) {
201 PathBuilder pb;
202 pb.moveTo(B.at0());
203 build_from_sbasis(pb, B, tol);
204 pb.finish();
205 return pb.peek().front();
206 }
208 //TODO: some of this logic should be lifted into svg-path
209 std::vector<Geom::Path>
210 path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol) {
211 Geom::PathBuilder pb;
212 if(B.size() == 0) return pb.peek();
213 Geom::Point start = B[0].at0();
214 pb.moveTo(start);
215 for(unsigned i = 0; ; i++) {
216 if(i+1 == B.size() || !are_near(B[i+1].at0(), B[i].at1(), tol)) {
217 //start of a new path
218 if(are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
219 //last line seg already there
220 goto no_add;
221 }
222 build_from_sbasis(pb, B[i], tol);
223 if(are_near(start, B[i].at1())) {
224 //it's closed
225 pb.closePath();
226 }
227 no_add:
228 if(i+1 >= B.size()) break;
229 start = B[i+1].at0();
230 pb.moveTo(start);
231 } else {
232 build_from_sbasis(pb, B[i], tol);
233 }
234 }
235 pb.finish();
236 return pb.peek();
237 }
239 };
241 /*
242 Local Variables:
243 mode:c++
244 c-file-style:"stroustrup"
245 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
246 indent-tabs-mode:nil
247 fill-column:99
248 End:
249 */
250 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :