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"
2
* Symmetric Power Basis - Bernstein Basis conversion routines
5
* Marco Cecchetti <mrcekets at gmail.com>
6
* Nathan Hurst <njh@mail.csse.monash.edu.au>
8
* Copyright 2007-2008 authors
10
* This library is free software; you can redistribute it and/or
11
* modify it either under the terms of the GNU Lesser General Public
12
* License version 2.1 as published by the Free Software Foundation
13
* (the "LGPL") or, at your option, under the terms of the Mozilla
14
* Public License Version 1.1 (the "MPL"). If you do not alter this
15
* notice, a recipient may use your version of this file under either
16
* the MPL or the LGPL.
18
* You should have received a copy of the LGPL along with this library
19
* in the file COPYING-LGPL-2.1; if not, write to the Free Software
20
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
* You should have received a copy of the MPL along with this library
22
* in the file COPYING-MPL-1.1
24
* The contents of this file are subject to the Mozilla Public License
25
* Version 1.1 (the "License"); you may not use this file except in
26
* compliance with the License. You may obtain a copy of the License at
27
* http://www.mozilla.org/MPL/
29
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30
* OF ANY KIND, either express or implied. See the LGPL or the MPL for
31
* the specific language governing rights and limitations.
35
#include <2geom/sbasis-to-bezier.h>
37
#include <2geom/choose.h>
38
#include <2geom/svg-path.h>
39
#include <2geom/exception.h>
15
41
#include <iostream>
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);
50
* Symmetric Power Basis - Bernstein Basis conversion routines
52
* some remark about precision:
53
* interval [0,1], subdivisions: 10^3
54
* - bezier_to_sbasis : up to degree ~72 precision is at least 10^-5
55
* up to degree ~87 precision is at least 10^-3
56
* - sbasis_to_bezier : up to order ~63 precision is at least 10^-15
57
* precision is at least 10^-14 even beyond order 200
59
* interval [-1,1], subdivisions: 10^3
60
* - bezier_to_sbasis : up to degree ~21 precision is at least 10^-5
61
* up to degree ~24 precision is at least 10^-3
62
* - sbasis_to_bezier : up to order ~11 precision is at least 10^-5
63
* up to order ~13 precision is at least 10^-3
65
* interval [-10,10], subdivisions: 10^3
66
* - bezier_to_sbasis : up to degree ~7 precision is at least 10^-5
67
* up to degree ~8 precision is at least 10^-3
68
* - sbasis_to_bezier : up to order ~3 precision is at least 10^-5
69
* up to order ~4 precision is at least 10^-3
72
* this implementation is based on the following article:
73
* J.Sanchez-Reyes - The Symmetric Analogue of the Polynomial Power Basis
77
double binomial(unsigned int n, unsigned int k)
79
return choose<double>(n, k);
83
int sgn(unsigned int j, unsigned int k)
86
// we are sure that j >= k
87
return ((j-k) & 1u) ? -1 : 1;
91
/** Changes the basis of p to be bernstein.
92
\param p the Symmetric basis polynomial
93
\returns the Bernstein basis polynomial
95
if the degree is even q is the order in the symmetrical power basis,
96
if the degree is odd q is the order + 1
97
n is always the polynomial degree, i. e. the Bezier order
98
sz is the number of bezier handles.
100
void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz)
107
if (sb[q-1][0] == sb[q-1][1])
121
q = (sz > 2*sb.size()-1) ? sb.size() : (sz+1)/2;
128
for (size_t k = 0; k < q; ++k)
130
for (size_t j = k; j < n-k; ++j) // j <= n-k-1
132
Tjk = binomial(n-2*k-1, j-k);
133
bz[j] += (Tjk * sb[k][0]);
134
bz[n-j] += (Tjk * sb[k][1]); // n-k <-> [k][1]
141
// the resulting coefficients are with respect to the scaled Bernstein
142
// basis so we need to divide them by (n, j) binomial coefficient
143
for (size_t j = 1; j < n; ++j)
145
bz[j] /= binomial(n, j);
151
/** Changes the basis of p to be Bernstein.
152
\param p the D2 Symmetric basis polynomial
153
\returns the D2 Bernstein basis polynomial
155
sz is always the polynomial degree, i. e. the Bezier order
157
void sbasis_to_bezier (std::vector<Point> & bz, D2<SBasis> const& sb, size_t sz)
161
sz = std::max(sb[X].size(), sb[Y].size())*2;
163
sbasis_to_bezier(bzx, sb[X], sz);
164
sbasis_to_bezier(bzy, sb[Y], sz);
165
assert(bzx.size() == bzy.size());
166
size_t n = (bzx.size() >= bzy.size()) ? bzx.size() : bzy.size();
168
bz.resize(n, Point(0,0));
169
for (size_t i = 0; i < bzx.size(); ++i)
173
for (size_t i = 0; i < bzy.size(); ++i)
180
/** Changes the basis of p to be sbasis.
181
\param p the Bernstein basis polynomial
182
\returns the Symmetric basis polynomial
184
if the degree is even q is the order in the symmetrical power basis,
185
if the degree is odd q is the order + 1
186
n is always the polynomial degree, i. e. the Bezier order
188
void bezier_to_sbasis (SBasis & sb, Bezier const& bz)
190
size_t n = bz.order();
191
size_t q = (n+1) / 2;
192
size_t even = (n & 1u) ? 0 : 1;
194
sb.resize(q + even, Linear(0, 0));
196
for (size_t k = 0; k < q; ++k)
198
for (size_t j = k; j < q; ++j)
200
Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
201
sb[j][0] += (Tjk * bz[k]);
202
sb[j][1] += (Tjk * bz[n-k]); // n-j <-> [j][1]
204
for (size_t j = k+1; j < q; ++j)
206
Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
207
sb[j][0] += (Tjk * bz[n-k]);
208
sb[j][1] += (Tjk * bz[k]); // n-j <-> [j][1]
213
for (size_t k = 0; k < q; ++k)
215
Tjk = sgn(q,k) * binomial(n, k);
216
sb[q][0] += (Tjk * (bz[k] + bz[n-k]));
218
sb[q][0] += (binomial(n, q) * bz[q]);
226
/** Changes the basis of d2 p to be sbasis.
227
\param p the d2 Bernstein basis polynomial
228
\returns the d2 Symmetric basis polynomial
230
if the degree is even q is the order in the symmetrical power basis,
231
if the degree is odd q is the order + 1
232
n is always the polynomial degree, i. e. the Bezier order
234
void bezier_to_sbasis (D2<SBasis> & sb, std::vector<Point> const& bz)
236
size_t n = bz.size() - 1;
237
size_t q = (n+1) / 2;
238
size_t even = (n & 1u) ? 0 : 1;
241
sb[X].resize(q + even, Linear(0, 0));
242
sb[Y].resize(q + even, Linear(0, 0));
244
for (size_t k = 0; k < q; ++k)
246
for (size_t j = k; j < q; ++j)
248
Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
249
sb[X][j][0] += (Tjk * bz[k][X]);
250
sb[X][j][1] += (Tjk * bz[n-k][X]);
251
sb[Y][j][0] += (Tjk * bz[k][Y]);
252
sb[Y][j][1] += (Tjk * bz[n-k][Y]);
254
for (size_t j = k+1; j < q; ++j)
256
Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
257
sb[X][j][0] += (Tjk * bz[n-k][X]);
258
sb[X][j][1] += (Tjk * bz[k][X]);
259
sb[Y][j][0] += (Tjk * bz[n-k][Y]);
260
sb[Y][j][1] += (Tjk * bz[k][Y]);
265
for (size_t k = 0; k < q; ++k)
267
Tjk = sgn(q,k) * binomial(n, k);
268
sb[X][q][0] += (Tjk * (bz[k][X] + bz[n-k][X]));
269
sb[Y][q][0] += (Tjk * (bz[k][Y] + bz[n-k][Y]));
271
sb[X][q][0] += (binomial(n, q) * bz[q][X]);
272
sb[X][q][1] = sb[X][q][0];
273
sb[Y][q][0] += (binomial(n, q) * bz[q][Y]);
274
sb[Y][q][1] = sb[Y][q][0];
276
sb[X][0][0] = bz[0][X];
277
sb[X][0][1] = bz[n][X];
278
sb[Y][0][0] = bz[0][Y];
279
sb[Y][0][1] = bz[n][Y];
283
} // end namespace Geom
136
288
* This version works by inverting a reasonable upper bound on the error term after subdividing the
137
289
* curve at $a$. We keep biting off pieces until there is no more curve left.
139
291
* Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$. A
140
292
* subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$. Let this be the desired
141
293
* tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$
171
323
pb.push_cubic(bez[1], bez[2], bez[3]);
173
325
// move to next piece of curve
174
326
if(a >= 1) break;
175
B = compose(B, Linear(a, 1));
327
B = compose(B, Linear(a, 1));
176
328
te = B.tail_error(k);
182
void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol) {
336
/** Make a path from a d2 sbasis.
337
\param p the d2 Symmetric basis polynomial
340
If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
342
void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
183
343
if (!B.isFinite()) {
184
throwException("assertion failed: B.isFinite()");
344
THROW_EXCEPTION("assertion failed: B.isFinite()");
186
346
if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
187
if(sbasis_size(B) <= 1) {
347
if( !only_cubicbeziers && (sbasis_size(B) <= 1) ) {
188
348
pb.lineTo(B.at1());
190
std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
350
std::vector<Geom::Point> bez;
351
sbasis_to_bezier(bez, B, 4);
191
352
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);
355
build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol, only_cubicbeziers);
356
build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol, only_cubicbeziers);
360
/** Make a path from a d2 sbasis.
361
\param p the d2 Symmetric basis polynomial
364
If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
200
path_from_sbasis(D2<SBasis> const &B, double tol) {
367
path_from_sbasis(D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
202
369
pb.moveTo(B.at0());
203
build_from_sbasis(pb, B, tol);
370
build_from_sbasis(pb, B, tol, only_cubicbeziers);
205
372
return pb.peek().front();
208
//TODO: some of this logic should be lifted into svg-path
375
/** Make a path from a d2 sbasis.
376
\param p the d2 Symmetric basis polynomial
379
If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
380
TODO: some of this logic should be lifted into svg-path
209
382
std::vector<Geom::Path>
210
path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol) {
383
path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol, bool only_cubicbeziers) {
211
384
Geom::PathBuilder pb;
212
385
if(B.size() == 0) return pb.peek();
213
386
Geom::Point start = B[0].at0();
added
removed
src/2geom/sbasis-to-bezier.cpp
