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
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
12
#include "sbasis-to-bezier.h"
16
#include "exception.h"
20
double W(unsigned n, unsigned j, unsigned k) {
22
if((n & 1) == 0 && j == q && k == q)
24
if(k > n-k) return W(n, n-j, n-k);
27
//assert(!(j >= n-k));
28
if(j >= n-k) return 0;
31
return choose<double>(n-2*k-1, j-k) /
35
// this produces a degree 2q bezier from a degree k sbasis
37
sbasis_to_bezier(SBasis const &B, unsigned q) {
40
/*if(B.back()[0] == B.back()[1]) {
45
Bezier result = Bezier(Bezier::Order(n-1));
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]);
62
// this produces a degree k sbasis from a degree 2q bezier
64
bezier_to_sbasis(Bezier const &B) {
65
unsigned n = B.size();
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]);
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;
88
result.resize(n, Geom::Point(0,0));
90
for(unsigned dim = 0; dim < 2; dim++) {
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]);
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]));
110
#if 0 // using old path
111
//std::vector<Geom::Point>
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
119
pb.start_subpath(Geom::Point(B[0][0][0], B[1][0][0]));
121
pb.push_line(Geom::Point(B[0][0][1], B[1][0][1]));
123
std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
125
pb.start_subpath(bez[0]);
127
pb.push_cubic(bez[1], bez[2], bez[3]);
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);
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.
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}$
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;
152
double A = std::sqrt(tol/te); // pow(te, 1./k)
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
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());
168
pb.start_subpath(bez[0]);
171
pb.push_cubic(bez[1], bez[2], bez[3]);
173
// move to next piece of curve
175
B = compose(B, Linear(a, 1));
176
te = B.tail_error(k);
182
void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol) {
184
throwException("assertion failed: B.isFinite()");
186
if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
187
if(sbasis_size(B) <= 1) {
190
std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
191
pb.curveTo(bez[1], bez[2], bez[3]);
194
build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol);
195
build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol);
200
path_from_sbasis(D2<SBasis> const &B, double tol) {
203
build_from_sbasis(pb, B, tol);
205
return pb.peek().front();
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();
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
222
build_from_sbasis(pb, B[i], tol);
223
if(are_near(start, B[i].at1())) {
228
if(i+1 >= B.size()) break;
229
start = B[i+1].at0();
232
build_from_sbasis(pb, B[i], tol);
244
c-file-style:"stroustrup"
245
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
250
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :