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

   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)
   5
   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)
  10
  11 */
  12 #include "sbasis-to-bezier.h"
  13 #include "choose.h"
  14 #include "svg-path.h"
  15 #include <iostream>
  16 #include "exception.h"
  17
  18 namespace Geom{
  19
  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 }
  34
  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 }
  57
  58 double mopi(int i) {
  59     return (i&1)?-1:1;
  60 }
  61
  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 }
  79
  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 */
 109
 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 }
 134
 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());
 149
 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]);
 172
 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 }
 179
 180 #endif
 181
 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 }
 198
 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 }
 207
 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 }
 238
 239 };
 240
 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 :