src/2geom/sbasis-to-bezier.cpp

   1 /*
   2  * Symmetric Power Basis - Bernstein Basis conversion routines
   3  *
   4  * Authors:
   5  *      Marco Cecchetti <mrcekets at gmail.com>
   6  *      Nathan Hurst <njh@mail.csse.monash.edu.au>
   7  *
   8  * Copyright 2007-2008  authors
   9  *
  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.
  17  *
  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
  23  *
  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/
  28  *
  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.
  32  */
  33
  34
  35 #include <2geom/sbasis-to-bezier.h>
  36 #include <2geom/d2.h>
  37 #include <2geom/choose.h>
  38 #include <2geom/svg-path.h>
  39 #include <2geom/exception.h>
  40
  41 #include <iostream>
  42
  43
  44
  45
  46 namespace Geom
  47 {
  48
  49 /*
  50  *  Symmetric Power Basis - Bernstein Basis conversion routines
  51  *
  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
  58  *
  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
  64  *
  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
  70  *
  71  *  references:
  72  *  this implementation is based on the following article:
  73  *  J.Sanchez-Reyes - The Symmetric Analogue of the Polynomial Power Basis
  74  */
  75
  76 inline
  77 double binomial(unsigned int n, unsigned int k)
  78 {
  79     return choose<double>(n, k);
  80 }
  81
  82 inline
  83 int sgn(unsigned int j, unsigned int k)
  84 {
  85     assert (j >= k);
  86     // we are sure that j >= k
  87     return ((j-k) &  1u) ? -1 : 1;
  88 }
  89
  90
  91 /** Changes the basis of p to be bernstein.
  92  \param p the Symmetric basis polynomial
  93  \returns the Bernstein basis polynomial
  94
  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.
  99 */
 100 void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz)
 101 {
 102     if (sb.size() == 0) {
 103         THROW_RANGEERROR("size of sb is too small");
 104     }
 105
 106     size_t q, n;
 107     bool even;
 108     if (sz == 0)
 109     {
 110         q = sb.size();
 111         if (sb[q-1][0] == sb[q-1][1])
 112         {
 113             even = true;
 114             --q;
 115             n = 2*q;
 116         }
 117         else
 118         {
 119             even = false;
 120             n = 2*q-1;
 121         }
 122     }
 123     else
 124     {
 125         q = (sz > 2*sb.size()-1) ?  sb.size() : (sz+1)/2;
 126         n = sz-1;
 127         even = false;
 128     }
 129     bz.clear();
 130     bz.resize(n+1);
 131     double Tjk;
 132     for (size_t k = 0; k < q; ++k)
 133     {
 134         for (size_t j = k; j < n-k; ++j) // j <= n-k-1
 135         {
 136             Tjk = binomial(n-2*k-1, j-k);
 137             bz[j] += (Tjk * sb[k][0]);
 138             bz[n-j] += (Tjk * sb[k][1]); // n-k <-> [k][1]
 139         }
 140     }
 141     if (even)
 142     {
 143         bz[q] += sb[q][0];
 144     }
 145     // the resulting coefficients are with respect to the scaled Bernstein
 146     // basis so we need to divide them by (n, j) binomial coefficient
 147     for (size_t j = 1; j < n; ++j)
 148     {
 149         bz[j] /= binomial(n, j);
 150     }
 151     bz[0] = sb[0][0];
 152     bz[n] = sb[0][1];
 153 }
 154
 155 /** Changes the basis of p to be Bernstein.
 156  \param p the D2 Symmetric basis polynomial
 157  \returns the D2 Bernstein basis polynomial
 158
 159  sz is always the polynomial degree, i. e. the Bezier order
 160 */
 161 void sbasis_to_bezier (std::vector<Point> & bz, D2<SBasis> const& sb, size_t sz)
 162 {
 163     Bezier bzx, bzy;
 164     if(sz == 0) {
 165         sz = std::max(sb[X].size(), sb[Y].size())*2;
 166     }
 167     sbasis_to_bezier(bzx, sb[X], sz);
 168     sbasis_to_bezier(bzy, sb[Y], sz);
 169     assert(bzx.size() == bzy.size());
 170     size_t n = (bzx.size() >= bzy.size()) ? bzx.size() : bzy.size();
 171
 172     bz.resize(n, Point(0,0));
 173     for (size_t i = 0; i < bzx.size(); ++i)
 174     {
 175         bz[i][X] = bzx[i];
 176     }
 177     for (size_t i = 0; i < bzy.size(); ++i)
 178     {
 179         bz[i][Y] = bzy[i];
 180     }
 181 }
 182
 183
 184 /** Changes the basis of p to be sbasis.
 185  \param p the Bernstein basis polynomial
 186  \returns the Symmetric basis polynomial
 187
 188  if the degree is even q is the order in the symmetrical power basis,
 189  if the degree is odd q is the order + 1
 190  n is always the polynomial degree, i. e. the Bezier order
 191 */
 192 void bezier_to_sbasis (SBasis & sb, Bezier const& bz)
 193 {
 194     size_t n = bz.order();
 195     size_t q = (n+1) / 2;
 196     size_t even = (n & 1u) ? 0 : 1;
 197     sb.clear();
 198     sb.resize(q + even, Linear(0, 0));
 199     double Tjk;
 200     for (size_t k = 0; k < q; ++k)
 201     {
 202         for (size_t j = k; j < q; ++j)
 203         {
 204             Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
 205             sb[j][0] += (Tjk * bz[k]);
 206             sb[j][1] += (Tjk * bz[n-k]); // n-j <-> [j][1]
 207         }
 208         for (size_t j = k+1; j < q; ++j)
 209         {
 210             Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
 211             sb[j][0] += (Tjk * bz[n-k]);
 212             sb[j][1] += (Tjk * bz[k]);   // n-j <-> [j][1]
 213         }
 214     }
 215     if (even)
 216     {
 217         for (size_t k = 0; k < q; ++k)
 218         {
 219             Tjk = sgn(q,k) * binomial(n, k);
 220             sb[q][0] += (Tjk * (bz[k] + bz[n-k]));
 221         }
 222         sb[q][0] += (binomial(n, q) * bz[q]);
 223         sb[q][1] = sb[q][0];
 224     }
 225     sb[0][0] = bz[0];
 226     sb[0][1] = bz[n];
 227 }
 228
 229
 230 /** Changes the basis of d2 p to be sbasis.
 231  \param p the d2 Bernstein basis polynomial
 232  \returns the d2 Symmetric basis polynomial
 233
 234  if the degree is even q is the order in the symmetrical power basis,
 235  if the degree is odd q is the order + 1
 236  n is always the polynomial degree, i. e. the Bezier order
 237 */
 238 void bezier_to_sbasis (D2<SBasis> & sb, std::vector<Point> const& bz)
 239 {
 240     size_t n = bz.size() - 1;
 241     size_t q = (n+1) / 2;
 242     size_t even = (n & 1u) ? 0 : 1;
 243     sb[X].clear();
 244     sb[Y].clear();
 245     sb[X].resize(q + even, Linear(0, 0));
 246     sb[Y].resize(q + even, Linear(0, 0));
 247     double Tjk;
 248     for (size_t k = 0; k < q; ++k)
 249     {
 250         for (size_t j = k; j < q; ++j)
 251         {
 252             Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
 253             sb[X][j][0] += (Tjk * bz[k][X]);
 254             sb[X][j][1] += (Tjk * bz[n-k][X]);
 255             sb[Y][j][0] += (Tjk * bz[k][Y]);
 256             sb[Y][j][1] += (Tjk * bz[n-k][Y]);
 257         }
 258         for (size_t j = k+1; j < q; ++j)
 259         {
 260             Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
 261             sb[X][j][0] += (Tjk * bz[n-k][X]);
 262             sb[X][j][1] += (Tjk * bz[k][X]);
 263             sb[Y][j][0] += (Tjk * bz[n-k][Y]);
 264             sb[Y][j][1] += (Tjk * bz[k][Y]);
 265         }
 266     }
 267     if (even)
 268     {
 269         for (size_t k = 0; k < q; ++k)
 270         {
 271             Tjk = sgn(q,k) * binomial(n, k);
 272             sb[X][q][0] += (Tjk * (bz[k][X] + bz[n-k][X]));
 273             sb[Y][q][0] += (Tjk * (bz[k][Y] + bz[n-k][Y]));
 274         }
 275         sb[X][q][0] += (binomial(n, q) * bz[q][X]);
 276         sb[X][q][1] = sb[X][q][0];
 277         sb[Y][q][0] += (binomial(n, q) * bz[q][Y]);
 278         sb[Y][q][1] = sb[Y][q][0];
 279     }
 280     sb[X][0][0] = bz[0][X];
 281     sb[X][0][1] = bz[n][X];
 282     sb[Y][0][0] = bz[0][Y];
 283     sb[Y][0][1] = bz[n][Y];
 284 }
 285
 286
 287 }  // end namespace Geom
 288
 289
 290 #if 0
 291 /*
 292 * This version works by inverting a reasonable upper bound on the error term after subdividing the
 293 * curve at $a$.  We keep biting off pieces until there is no more curve left.
 294 *
 295 * Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$.  A
 296 * subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$.  Let this be the desired
 297 * tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$
 298 */
 299 void
 300 subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, double tol, bool initial) {
 301     const unsigned k = 2; // cubic bezier
 302     double te = B.tail_error(k);
 303     assert(B[0].IS_FINITE());
 304     assert(B[1].IS_FINITE());
 305
 306     //std::cout << "tol = " << tol << std::endl;
 307     while(1) {
 308         double A = std::sqrt(tol/te); // pow(te, 1./k)
 309         double a = A;
 310         if(A < 1) {
 311             A = std::min(A, 0.25);
 312             a = 0.5 - std::sqrt(0.25 - A); // quadratic formula
 313             if(a > 1) a = 1; // clamp to the end of the segment
 314         } else
 315             a = 1;
 316         assert(a > 0);
 317         //std::cout << "te = " << te << std::endl;
 318         //std::cout << "A = " << A << "; a=" << a << std::endl;
 319         D2<SBasis> Bs = compose(B, Linear(0, a));
 320         assert(Bs.tail_error(k));
 321         std::vector<Geom::Point> bez = sbasis_to_bezier(Bs, 2);
 322         reverse(bez.begin(), bez.end());
 323         if (initial) {
 324           pb.start_subpath(bez[0]);
 325           initial = false;
 326         }
 327         pb.push_cubic(bez[1], bez[2], bez[3]);
 328
 329 // move to next piece of curve
 330         if(a >= 1) break;
 331         B = compose(B, Linear(a, 1));
 332         te = B.tail_error(k);
 333     }
 334 }
 335
 336 #endif
 337
 338 namespace Geom{
 339
 340 /** Make a path from a d2 sbasis.
 341  \param p the d2 Symmetric basis polynomial
 342  \returns a Path
 343
 344   If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
 345 */
 346 void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
 347     if (!B.isFinite()) {
 348         THROW_EXCEPTION("assertion failed: B.isFinite()");
 349     }
 350     if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
 351         if( !only_cubicbeziers && (sbasis_size(B) <= 1) ) {
 352             pb.lineTo(B.at1());
 353         } else {
 354             std::vector<Geom::Point> bez;
 355             sbasis_to_bezier(bez, B, 4);
 356             pb.curveTo(bez[1], bez[2], bez[3]);
 357         }
 358     } else {
 359         build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol, only_cubicbeziers);
 360         build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol, only_cubicbeziers);
 361     }
 362 }
 363
 364 /** Make a path from a d2 sbasis.
 365  \param p the d2 Symmetric basis polynomial
 366  \returns a Path
 367
 368   If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
 369 */
 370 Path
 371 path_from_sbasis(D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
 372     PathBuilder pb;
 373     pb.moveTo(B.at0());
 374     build_from_sbasis(pb, B, tol, only_cubicbeziers);
 375     pb.finish();
 376     return pb.peek().front();
 377 }
 378
 379 /** Make a path from a d2 sbasis.
 380  \param p the d2 Symmetric basis polynomial
 381  \returns a Path
 382
 383   If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
 384  TODO: some of this logic should be lifted into svg-path
 385 */
 386 std::vector<Geom::Path>
 387 path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol, bool only_cubicbeziers) {
 388     Geom::PathBuilder pb;
 389     if(B.size() == 0) return pb.peek();
 390     Geom::Point start = B[0].at0();
 391     pb.moveTo(start);
 392     for(unsigned i = 0; ; i++) {
 393         if(i+1 == B.size() || !are_near(B[i+1].at0(), B[i].at1(), tol)) {
 394             //start of a new path
 395             if(are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
 396                 pb.closePath();
 397                 //last line seg already there (because of .closePath())
 398                 goto no_add;
 399             }
 400             build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
 401             if(are_near(start, B[i].at1())) {
 402                 //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.
 403                 pb.closePath();
 404             }
 405           no_add:
 406             if(i+1 >= B.size()) break;
 407             start = B[i+1].at0();
 408             pb.moveTo(start);
 409         } else {
 410             build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
 411         }
 412     }
 413     pb.finish();
 414     return pb.peek();
 415 }
 416
 417 }
 418
 419 /*
 420   Local Variables:
 421   mode:c++
 422   c-file-style:"stroustrup"
 423   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 424   indent-tabs-mode:nil
 425   fill-column:99
 426   End:
 427 */
 428 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :