Adding upstream version 1.0.2.upstream/1.0.2upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 570 insertions, 0 deletions

diff --git a/src/2geom/sbasis-to-bezier.cpp b/src/2geom/sbasis-to-bezier.cpp
new file mode 100644
index 0000000..010cf7c
--- /dev/null
+++ b/src/2geom/sbasis-to-bezier.cpp
@@ -0,0 +1,570 @@
+/*
+ * Symmetric Power Basis - Bernstein Basis conversion routines
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *      Nathan Hurst <njh@mail.csse.monash.edu.au>
+ *
+ * Copyright 2007-2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/d2.h>
+#include <2geom/choose.h>
+#include <2geom/path-sink.h>
+#include <2geom/exception.h>
+#include <2geom/convex-hull.h>
+
+#include <iostream>
+
+
+
+
+namespace Geom
+{
+
+/*
+ *  Symmetric Power Basis - Bernstein Basis conversion routines
+ *
+ *  some remark about precision:
+ *  interval [0,1], subdivisions: 10^3
+ *  - bezier_to_sbasis : up to degree ~72 precision is at least 10^-5
+ *                       up to degree ~87 precision is at least 10^-3
+ *  - sbasis_to_bezier : up to order ~63 precision is at least 10^-15
+ *                       precision is at least 10^-14 even beyond order 200
+ *
+ *  interval [-1,1], subdivisions: 10^3
+ *  - bezier_to_sbasis : up to degree ~21 precision is at least 10^-5
+ *                       up to degree ~24 precision is at least 10^-3
+ *  - sbasis_to_bezier : up to order ~11 precision is at least 10^-5
+ *                       up to order ~13 precision is at least 10^-3
+ *
+ *  interval [-10,10], subdivisions: 10^3
+ *  - bezier_to_sbasis : up to degree ~7 precision is at least 10^-5
+ *                       up to degree ~8 precision is at least 10^-3
+ *  - sbasis_to_bezier : up to order ~3 precision is at least 10^-5
+ *                       up to order ~4 precision is at least 10^-3
+ *
+ *  references:
+ *  this implementation is based on the following article:
+ *  J.Sanchez-Reyes - The Symmetric Analogue of the Polynomial Power Basis
+ */
+
+inline
+double binomial(unsigned int n, unsigned int k)
+{
+    return choose<double>(n, k);
+}
+
+inline
+int sgn(unsigned int j, unsigned int k)
+{
+    assert (j >= k);
+    // we are sure that j >= k
+    return ((j-k) &  1u) ? -1 : 1;
+}
+
+
+/** Changes the basis of p to be bernstein.
+ \param p the Symmetric basis polynomial
+ \returns the Bernstein basis polynomial
+
+ if the degree is even q is the order in the symmetrical power basis,
+ if the degree is odd q is the order + 1
+ n is always the polynomial degree, i. e. the Bezier order
+ sz is the number of bezier handles.
+*/
+void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz)
+{
+    assert(sb.size() > 0);
+
+    size_t q, n;
+    bool even;
+    if (sz == 0)
+    {
+        q = sb.size();
+        if (sb[q-1][0] == sb[q-1][1])
+        {
+            even = true;
+            --q;
+            n = 2*q;
+        }
+        else
+        {
+            even = false;
+            n = 2*q-1;
+        }
+    }
+    else
+    {
+        q = (sz > 2*sb.size()-1) ?  sb.size() : (sz+1)/2;
+        n = sz-1;
+        even = false;
+    }
+    bz.clear();
+    bz.resize(n+1);
+    double Tjk;
+    for (size_t k = 0; k < q; ++k)
+    {
+        for (size_t j = k; j < n-k; ++j) // j <= n-k-1
+        {
+            Tjk = binomial(n-2*k-1, j-k);
+            bz[j] += (Tjk * sb[k][0]);
+            bz[n-j] += (Tjk * sb[k][1]); // n-k <-> [k][1]
+        }
+    }
+    if (even)
+    {
+        bz[q] += sb[q][0];
+    }
+    // the resulting coefficients are with respect to the scaled Bernstein
+    // basis so we need to divide them by (n, j) binomial coefficient
+    for (size_t j = 1; j < n; ++j)
+    {
+        bz[j] /= binomial(n, j);
+    }
+    bz[0] = sb[0][0];
+    bz[n] = sb[0][1];
+}
+
+void sbasis_to_bezier(D2<Bezier> &bz, D2<SBasis> const &sb, size_t sz)
+{
+    if (sz == 0) {
+        sz = std::max(sb[X].size(), sb[Y].size())*2;
+    }
+    sbasis_to_bezier(bz[X], sb[X], sz);
+    sbasis_to_bezier(bz[Y], sb[Y], sz);
+}
+
+/** Changes the basis of p to be Bernstein.
+ \param p the D2 Symmetric basis polynomial
+ \returns the D2 Bernstein basis polynomial
+
+ sz is always the polynomial degree, i. e. the Bezier order
+*/
+void sbasis_to_bezier (std::vector<Point> & bz, D2<SBasis> const& sb, size_t sz)
+{
+    D2<Bezier> bez;
+    sbasis_to_bezier(bez, sb, sz);
+    bz = bezier_points(bez);
+}
+
+/** Changes the basis of p to be Bernstein.
+ \param p the D2 Symmetric basis polynomial
+ \returns the D2 Bernstein basis cubic polynomial
+
+Bezier is always cubic.
+For general asymmetric case, fit the SBasis function value at midpoint
+For parallel, symmetric case, find the point of closest approach to the midpoint
+For parallel, anti-symmetric case, fit the SBasis slope at midpoint
+*/
+void sbasis_to_cubic_bezier (std::vector<Point> & bz, D2<SBasis> const& sb)
+{
+    double delx[2], dely[2];
+    double xprime[2], yprime[2];
+    double midx = 0;
+    double midy = 0;
+    double midx_0, midy_0;
+    double numer[2], numer_0[2];
+    double denom;
+    double div;
+
+    if ((sb[X].size() == 0) || (sb[Y].size() == 0)) {
+        THROW_RANGEERROR("size of sb is too small");
+    }
+
+    sbasis_to_bezier(bz, sb, 4);  // zeroth-order estimate
+    if ((sb[X].size() < 3) && (sb[Y].size() < 3))
+        return;  // cubic bezier estimate is exact
+    Geom::ConvexHull bezhull(bz);
+
+//  calculate first derivatives of x and y wrt t
+
+    for (int i = 0; i < 2; ++i) {
+        xprime[i] = sb[X][0][1] - sb[X][0][0];
+        yprime[i] = sb[Y][0][1] - sb[Y][0][0];
+    }
+    if (sb[X].size() > 1) {
+        xprime[0] += sb[X][1][0];
+        xprime[1] -= sb[X][1][1];
+    }
+    if (sb[Y].size() > 1) {
+        yprime[0] += sb[Y][1][0];
+        yprime[1] -= sb[Y][1][1];
+    }
+
+//  calculate midpoint at t = 0.5
+
+    div = 2;
+    for (size_t i = 0; i < sb[X].size(); ++i) {
+        midx += (sb[X][i][0] + sb[X][i][1])/div;
+        div *= 4;
+    }
+
+    div = 2;
+    for (size_t i = 0; i < sb[Y].size(); ++i) {
+        midy += (sb[Y][i][0] + sb[Y][i][1])/div;
+        div *= 4;
+    }
+
+//  is midpoint in hull: if not, the solution will be ill-conditioned, LP Bug 1428683
+
+    if (!bezhull.contains(Geom::Point(midx, midy)))
+        return;
+
+//  calculate Bezier control arms
+
+    midx = 8*midx - 4*bz[0][X] - 4*bz[3][X];  // re-define relative to center
+    midy = 8*midy - 4*bz[0][Y] - 4*bz[3][Y];
+    midx_0 = sb[X][1][0] + sb[X][1][1];       // zeroth order estimate
+    midy_0 = sb[Y][1][0] + sb[Y][1][1];
+
+    if ((std::abs(xprime[0]) < EPSILON) && (std::abs(yprime[0]) < EPSILON)
+    && ((std::abs(xprime[1]) > EPSILON) || (std::abs(yprime[1]) > EPSILON)))  { // degenerate handle at 0 : use distance of closest approach
+        numer[0] = midx*xprime[1] + midy*yprime[1];
+        denom = 3.0*(xprime[1]*xprime[1] + yprime[1]*yprime[1]);
+        delx[0] = 0;
+        dely[0] = 0;
+        delx[1] = -xprime[1]*numer[0]/denom;
+        dely[1] = -yprime[1]*numer[0]/denom;
+    } else if ((std::abs(xprime[1]) < EPSILON) && (std::abs(yprime[1]) < EPSILON)
+           && ((std::abs(xprime[0]) > EPSILON) || (std::abs(yprime[0]) > EPSILON)))  { // degenerate handle at 1 : ditto
+        numer[1] = midx*xprime[0] + midy*yprime[0];
+        denom = 3.0*(xprime[0]*xprime[0] + yprime[0]*yprime[0]);
+        delx[0] = xprime[0]*numer[1]/denom;
+        dely[0] = yprime[0]*numer[1]/denom;
+        delx[1] = 0;
+        dely[1] = 0;
+    } else if  (std::abs(xprime[1]*yprime[0] - yprime[1]*xprime[0]) >  // general case : fit mid fxn value
+        0.002 * std::abs(xprime[1]*xprime[0] + yprime[1]*yprime[0])) { // approx. 0.1 degree of angle
+        double test1 = (bz[1][Y] - bz[0][Y])*(bz[3][X] - bz[0][X]) - (bz[1][X] - bz[0][X])*(bz[3][Y] - bz[0][Y]);
+        double test2 = (bz[2][Y] - bz[0][Y])*(bz[3][X] - bz[0][X]) - (bz[2][X] - bz[0][X])*(bz[3][Y] - bz[0][Y]);
+        if (test1*test2 < 0) // reject anti-symmetric case, LP Bug 1428267 & Bug 1428683
+            return;
+        denom = 3.0*(xprime[1]*yprime[0] - yprime[1]*xprime[0]);
+        for (int i = 0; i < 2; ++i) {
+            numer_0[i] = xprime[1 - i]*midy_0 - yprime[1 - i]*midx_0;
+            numer[i] = xprime[1 - i]*midy - yprime[1 - i]*midx;
+            delx[i] = xprime[i]*numer[i]/denom;
+            dely[i] = yprime[i]*numer[i]/denom;
+            if (numer_0[i]*numer[i] < 0)  // check for reversal of direction, LP Bug 1544680
+                return;
+        }
+        if (std::abs((numer[0] - numer_0[0])*numer_0[1]) > 10.0*std::abs((numer[1] - numer_0[1])*numer_0[0]) // check for asymmetry
+        ||  std::abs((numer[1] - numer_0[1])*numer_0[0]) > 10.0*std::abs((numer[0] - numer_0[0])*numer_0[1]))
+            return;
+    } else if ((xprime[0]*xprime[1] < 0) || (yprime[0]*yprime[1] < 0)) { // symmetric case : use distance of closest approach
+        numer[0] = midx*xprime[0] + midy*yprime[0];
+        denom = 6.0*(xprime[0]*xprime[0] + yprime[0]*yprime[0]);
+        delx[0] = xprime[0]*numer[0]/denom;
+        dely[0] = yprime[0]*numer[0]/denom;
+        delx[1] = -delx[0];
+        dely[1] = -dely[0];
+    } else {                                    // anti-symmetric case : fit mid slope
+                                                // calculate slope at t = 0.5
+        midx = 0;
+        div = 1;
+        for (size_t i = 0; i < sb[X].size(); ++i) {
+            midx += (sb[X][i][1] - sb[X][i][0])/div;
+            div *= 4;
+        }
+        midy = 0;
+        div = 1;
+        for (size_t i = 0; i < sb[Y].size(); ++i) {
+            midy += (sb[Y][i][1] - sb[Y][i][0])/div;
+            div *= 4;
+        }
+        if (midx*yprime[0] != midy*xprime[0]) {
+            denom = midx*yprime[0] - midy*xprime[0];
+            numer[0] = midx*(bz[3][Y] - bz[0][Y]) - midy*(bz[3][X] - bz[0][X]);
+            for (int i = 0; i < 2; ++i) {
+                delx[i] = xprime[0]*numer[0]/denom;
+                dely[i] = yprime[0]*numer[0]/denom;
+            }
+        } else {                                // linear case
+            for (int i = 0; i < 2; ++i) {
+                delx[i] = (bz[3][X] - bz[0][X])/3;
+                dely[i] = (bz[3][Y] - bz[0][Y])/3;
+            }
+        }
+    }
+    bz[1][X] = bz[0][X] + delx[0];
+    bz[1][Y] = bz[0][Y] + dely[0];
+    bz[2][X] = bz[3][X] - delx[1];
+    bz[2][Y] = bz[3][Y] - dely[1];
+}
+
+/** Changes the basis of p to be sbasis.
+ \param p the Bernstein basis polynomial
+ \returns the Symmetric basis polynomial
+
+ if the degree is even q is the order in the symmetrical power basis,
+ if the degree is odd q is the order + 1
+ n is always the polynomial degree, i. e. the Bezier order
+*/
+void bezier_to_sbasis (SBasis & sb, Bezier const& bz)
+{
+    size_t n = bz.order();
+    size_t q = (n+1) / 2;
+    size_t even = (n & 1u) ? 0 : 1;
+    sb.clear();
+    sb.resize(q + even, Linear(0, 0));
+    double Tjk;
+    for (size_t k = 0; k < q; ++k)
+    {
+        for (size_t j = k; j < q; ++j)
+        {
+            Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
+            sb[j][0] += (Tjk * bz[k]);
+            sb[j][1] += (Tjk * bz[n-k]); // n-j <-> [j][1]
+        }
+        for (size_t j = k+1; j < q; ++j)
+        {
+            Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
+            sb[j][0] += (Tjk * bz[n-k]);
+            sb[j][1] += (Tjk * bz[k]);   // n-j <-> [j][1]
+        }
+    }
+    if (even)
+    {
+        for (size_t k = 0; k < q; ++k)
+        {
+            Tjk = sgn(q,k) * binomial(n, k);
+            sb[q][0] += (Tjk * (bz[k] + bz[n-k]));
+        }
+        sb[q][0] += (binomial(n, q) * bz[q]);
+        sb[q][1] = sb[q][0];
+    }
+    sb[0][0] = bz[0];
+    sb[0][1] = bz[n];
+}
+
+
+/** Changes the basis of d2 p to be sbasis.
+ \param p the d2 Bernstein basis polynomial
+ \returns the d2 Symmetric basis polynomial
+
+ if the degree is even q is the order in the symmetrical power basis,
+ if the degree is odd q is the order + 1
+ n is always the polynomial degree, i. e. the Bezier order
+*/
+void bezier_to_sbasis (D2<SBasis> & sb, std::vector<Point> const& bz)
+{
+    size_t n = bz.size() - 1;
+    size_t q = (n+1) / 2;
+    size_t even = (n & 1u) ? 0 : 1;
+    sb[X].clear();
+    sb[Y].clear();
+    sb[X].resize(q + even, Linear(0, 0));
+    sb[Y].resize(q + even, Linear(0, 0));
+    double Tjk;
+    for (size_t k = 0; k < q; ++k)
+    {
+        for (size_t j = k; j < q; ++j)
+        {
+            Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
+            sb[X][j][0] += (Tjk * bz[k][X]);
+            sb[X][j][1] += (Tjk * bz[n-k][X]);
+            sb[Y][j][0] += (Tjk * bz[k][Y]);
+            sb[Y][j][1] += (Tjk * bz[n-k][Y]);
+        }
+        for (size_t j = k+1; j < q; ++j)
+        {
+            Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
+            sb[X][j][0] += (Tjk * bz[n-k][X]);
+            sb[X][j][1] += (Tjk * bz[k][X]);
+            sb[Y][j][0] += (Tjk * bz[n-k][Y]);
+            sb[Y][j][1] += (Tjk * bz[k][Y]);
+        }
+    }
+    if (even)
+    {
+        for (size_t k = 0; k < q; ++k)
+        {
+            Tjk = sgn(q,k) * binomial(n, k);
+            sb[X][q][0] += (Tjk * (bz[k][X] + bz[n-k][X]));
+            sb[Y][q][0] += (Tjk * (bz[k][Y] + bz[n-k][Y]));
+        }
+        sb[X][q][0] += (binomial(n, q) * bz[q][X]);
+        sb[X][q][1] = sb[X][q][0];
+        sb[Y][q][0] += (binomial(n, q) * bz[q][Y]);
+        sb[Y][q][1] = sb[Y][q][0];
+    }
+    sb[X][0][0] = bz[0][X];
+    sb[X][0][1] = bz[n][X];
+    sb[Y][0][0] = bz[0][Y];
+    sb[Y][0][1] = bz[n][Y];
+}
+
+
+}  // end namespace Geom
+
+
+#if 0 
+/*
+* This version works by inverting a reasonable upper bound on the error term after subdividing the
+* curve at $a$.  We keep biting off pieces until there is no more curve left.
+*
+* Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$.  A
+* subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$.  Let this be the desired
+* tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$
+*/
+void
+subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, double tol, bool initial) {
+    const unsigned k = 2; // cubic bezier
+    double te = B.tail_error(k);
+    assert(B[0].std::isfinite());
+    assert(B[1].std::isfinite());
+
+    //std::cout << "tol = " << tol << std::endl;
+    while(1) {
+        double A = std::sqrt(tol/te); // pow(te, 1./k)
+        double a = A;
+        if(A < 1) {
+            A = std::min(A, 0.25);
+            a = 0.5 - std::sqrt(0.25 - A); // quadratic formula
+            if(a > 1) a = 1; // clamp to the end of the segment
+        } else
+            a = 1;
+        assert(a > 0);
+        //std::cout << "te = " << te << std::endl;
+        //std::cout << "A = " << A << "; a=" << a << std::endl;
+        D2<SBasis> Bs = compose(B, Linear(0, a));
+        assert(Bs.tail_error(k));
+        std::vector<Geom::Point> bez = sbasis_to_bezier(Bs, 2);
+        reverse(bez.begin(), bez.end());
+        if (initial) {
+          pb.start_subpath(bez[0]);
+          initial = false;
+        }
+        pb.push_cubic(bez[1], bez[2], bez[3]);
+
+// move to next piece of curve
+        if(a >= 1) break;
+        B = compose(B, Linear(a, 1));
+        te = B.tail_error(k);
+    }
+}
+
+#endif
+
+namespace Geom{
+
+/** Make a path from a d2 sbasis.
+ \param p the d2 Symmetric basis polynomial
+ \returns a Path
+
+  If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
+*/
+void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
+    if (!B.isFinite()) {
+        THROW_EXCEPTION("assertion failed: B.isFinite()");
+    }
+    if(tail_error(B, 3) < tol || sbasis_size(B) == 2) { // nearly cubic enough
+        if( !only_cubicbeziers && (sbasis_size(B) <= 1) ) {
+            pb.lineTo(B.at1());
+        } else {
+            std::vector<Geom::Point> bez;
+//            sbasis_to_bezier(bez, B, 4);
+            sbasis_to_cubic_bezier(bez, B);
+            pb.curveTo(bez[1], bez[2], bez[3]);
+        }
+    } else {
+        build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol, only_cubicbeziers);
+        build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol, only_cubicbeziers);
+    }
+}
+
+/** Make a path from a d2 sbasis.
+ \param p the d2 Symmetric basis polynomial
+ \returns a Path
+
+  If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
+*/
+Path
+path_from_sbasis(D2<SBasis> const &B, double tol, bool only_cubicbeziers) {
+    PathBuilder pb;
+    pb.moveTo(B.at0());
+    build_from_sbasis(pb, B, tol, only_cubicbeziers);
+    pb.flush();
+    return pb.peek().front();
+}
+
+/** Make a path from a d2 sbasis.
+ \param p the d2 Symmetric basis polynomial
+ \returns a Path
+
+  If only_cubicbeziers is true, the resulting path may only contain CubicBezier curves.
+ TODO: some of this logic should be lifted into svg-path
+*/
+PathVector
+path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol, bool only_cubicbeziers) {
+    Geom::PathBuilder pb;
+    if(B.size() == 0) return pb.peek();
+    Geom::Point start = B[0].at0();
+    pb.moveTo(start);
+    for(unsigned i = 0; ; i++) {
+        if ( (i+1 == B.size()) 
+             || !are_near(B[i+1].at0(), B[i].at1(), tol) )
+        {
+            //start of a new path
+            if (are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
+                pb.closePath();
+                //last line seg already there (because of .closePath())
+                goto no_add;
+            }
+            build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
+            if (are_near(start, B[i].at1())) {
+                //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.
+                pb.closePath();
+            }
+          no_add:
+            if (i+1 >= B.size()) {
+                break;
+            }
+            start = B[i+1].at0();
+            pb.moveTo(start);
+        } else {
+            build_from_sbasis(pb, B[i], tol, only_cubicbeziers);
+        }
+    }
+    pb.flush();
+    return pb.peek();
+}
+
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :