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
  17 namespace Geom{
  18
  19 double W(unsigned n, unsigned j, unsigned k) {
  20     unsigned q = (n+1)/2;
  21     if((n & 1) == 0 && j == q && k == q)
  22         return 1;
  23     if(k > n-k) return W(n, n-j, n-k);
  24     assert((k <= q));
  25     if(k >= q) return 0;
  26     //assert(!(j >= n-k));
  27     if(j >= n-k) return 0;
  28     //assert(!(j < k));
  29     if(j < k) return 0;
  30     return choose<double>(n-2*k-1, j-k) /
  31         choose<double>(n,j);
  32 }
  33
  34 // this produces a degree 2q bezier from a degree k sbasis
  35 Bezier
  36 sbasis_to_bezier(SBasis const &B, unsigned q) {
  37     if(q == 0) {
  38         q = B.size();
  39         /*if(B.back()[0] == B.back()[1]) {
  40             n--;
  41             }*/
  42     }
  43     unsigned n = q*2;
  44     Bezier result = Bezier(Bezier::Order(n-1));
  45     if(q > B.size())
  46         q = B.size();
  47     n--;
  48     for(unsigned k = 0; k < q; k++) {
  49         for(unsigned j = 0; j <= n-k; j++) {
  50             result[j] += (W(n, j, k)*B[k][0] +
  51                           W(n, n-j, k)*B[k][1]);
  52         }
  53     }
  54     return result;
  55 }
  56
  57 double mopi(int i) {
  58     return (i&1)?-1:1;
  59 }
  60
  61 // this produces a degree k sbasis from a degree 2q bezier
  62 SBasis
  63 bezier_to_sbasis(Bezier const &B) {
  64     unsigned n = B.size();
  65     unsigned q = (n+1)/2;
  66     SBasis result;
  67     result.resize(q+1);
  68     for(unsigned k = 0; k < q; k++) {
  69         result[k][0] = result[k][1] = 0;
  70         for(unsigned j = 0; j <= n-k; j++) {
  71             result[k][0] += mopi(int(j)-int(k))*W(n, j, k)*B[j];
  72             result[k][1] += mopi(int(j)-int(k))*W(n, j, k)*B[j];
  73             //W(n, n-j, k)*B[k][1]);
  74         }
  75     }
  76     return result;
  77 }
  78
  79 // this produces a 2q point bezier from a degree q sbasis
  80 std::vector<Geom::Point>
  81 sbasis_to_bezier(D2<SBasis> const &B, unsigned qq) {
  82     std::vector<Geom::Point> result;
  83     if(qq == 0) {
  84         qq = sbasis_size(B);
  85     }
  86     unsigned n = qq * 2;
  87     result.resize(n, Geom::Point(0,0));
  88     n--;
  89     for(unsigned dim = 0; dim < 2; dim++) {
  90         unsigned q = qq;
  91         if(q > B[dim].size())
  92             q = B[dim].size();
  93         for(unsigned k = 0; k < q; k++) {
  94             for(unsigned j = 0; j <= n-k; j++) {
  95                 result[j][dim] += (W(n, j, k)*B[dim][k][0] +
  96                              W(n, n-j, k)*B[dim][k][1]);
  97                 }
  98         }
  99     }
 100     return result;
 101 }
 102 /*
 103 template <unsigned order>
 104 D2<Bezier<order> > sbasis_to_bezier(D2<SBasis> const &B) {
 105     return D2<Bezier<order> >(sbasis_to_bezier<order>(B[0]), sbasis_to_bezier<order>(B[1]));
 106 }
 107 */
 108
 109 #if 0 // using old path
 110 //std::vector<Geom::Point>
 111 // mutating
 112 void
 113 subpath_from_sbasis(Geom::OldPathSetBuilder &pb, D2<SBasis> const &B, double tol, bool initial) {
 114     assert(B.is_finite());
 115     if(B.tail_error(2) < tol || B.size() == 2) { // nearly cubic enough
 116         if(B.size() == 1) {
 117             if (initial) {
 118                 pb.start_subpath(Geom::Point(B[0][0][0], B[1][0][0]));
 119             }
 120             pb.push_line(Geom::Point(B[0][0][1], B[1][0][1]));
 121         } else {
 122             std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
 123             if (initial) {
 124                 pb.start_subpath(bez[0]);
 125             }
 126             pb.push_cubic(bez[1], bez[2], bez[3]);
 127         }
 128     } else {
 129         subpath_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol, initial);
 130         subpath_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol, false);
 131     }
 132 }
 133
 134 /*
 135 * This version works by inverting a reasonable upper bound on the error term after subdividing the
 136 * curve at $a$.  We keep biting off pieces until there is no more curve left.
 137 *
 138 * Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$.  A
 139 * subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$.  Let this be the desired
 140 * tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$
 141 */
 142 void
 143 subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, double tol, bool initial) {
 144     const unsigned k = 2; // cubic bezier
 145     double te = B.tail_error(k);
 146     assert(B[0].is_finite());
 147     assert(B[1].is_finite());
 148
 149     //std::cout << "tol = " << tol << std::endl;
 150     while(1) {
 151         double A = std::sqrt(tol/te); // pow(te, 1./k)
 152         double a = A;
 153         if(A < 1) {
 154             A = std::min(A, 0.25);
 155             a = 0.5 - std::sqrt(0.25 - A); // quadratic formula
 156             if(a > 1) a = 1; // clamp to the end of the segment
 157         } else
 158             a = 1;
 159         assert(a > 0);
 160         //std::cout << "te = " << te << std::endl;
 161         //std::cout << "A = " << A << "; a=" << a << std::endl;
 162         D2<SBasis> Bs = compose(B, Linear(0, a));
 163         assert(Bs.tail_error(k));
 164         std::vector<Geom::Point> bez = sbasis_to_bezier(Bs, 2);
 165         reverse(bez.begin(), bez.end());
 166         if (initial) {
 167           pb.start_subpath(bez[0]);
 168           initial = false;
 169         }
 170         pb.push_cubic(bez[1], bez[2], bez[3]);
 171
 172 // move to next piece of curve
 173         if(a >= 1) break;
 174         B = compose(B, Linear(a, 1));
 175         te = B.tail_error(k);
 176     }
 177 }
 178
 179 #endif
 180
 181 void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol) {
 182     assert(B.isFinite());
 183     if(tail_error(B, 2) < tol || sbasis_size(B) == 2) { // nearly cubic enough
 184         if(sbasis_size(B) <= 1) {
 185             pb.lineTo(B.at1());
 186         } else {
 187             std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
 188             pb.curveTo(bez[1], bez[2], bez[3]);
 189         }
 190     } else {
 191         build_from_sbasis(pb, compose(B, Linear(0, 0.5)), tol);
 192         build_from_sbasis(pb, compose(B, Linear(0.5, 1)), tol);
 193     }
 194 }
 195
 196 Path
 197 path_from_sbasis(D2<SBasis> const &B, double tol) {
 198     PathBuilder pb;
 199     pb.moveTo(B.at0());
 200     build_from_sbasis(pb, B, tol);
 201     pb.finish();
 202     return pb.peek().front();
 203 }
 204
 205 //TODO: some of this logic should be lifted into svg-path
 206 std::vector<Geom::Path>
 207 path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double tol) {
 208     Geom::PathBuilder pb;
 209     if(B.size() == 0) return pb.peek();
 210     Geom::Point start = B[0].at0();
 211     pb.moveTo(start);
 212     for(unsigned i = 0; ; i++) {
 213         if(i+1 == B.size() || !are_near(B[i+1].at0(), B[i].at1(), tol)) {
 214             //start of a new path
 215             if(are_near(start, B[i].at1()) && sbasis_size(B[i]) <= 1) {
 216                 //last line seg already there
 217                 goto no_add;
 218             }
 219             build_from_sbasis(pb, B[i], tol);
 220             if(are_near(start, B[i].at1())) {
 221                 //it's closed
 222                 pb.closePath();
 223             }
 224           no_add:
 225             if(i+1 >= B.size()) break;
 226             start = B[i+1].at0();
 227             pb.moveTo(start);
 228         } else {
 229             build_from_sbasis(pb, B[i], tol);
 230         }
 231     }
 232     pb.finish();
 233     return pb.peek();
 234 }
 235
 236 };
 237
 238 /*
 239   Local Variables:
 240   mode:c++
 241   c-file-style:"stroustrup"
 242   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 243   indent-tabs-mode:nil
 244   fill-column:99
 245   End:
 246 */
 247 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :