~ubuntu-branches/debian/squeeze/inkscape/squeeze

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Thomas Viehmann
  • Date: 2008-09-09 23:29:02 UTC
  • mfrom: (1.1.7 upstream)
  • Revision ID: james.westby@ubuntu.com-20080909232902-c50iujhk1w79u8e7
Tags: 0.46-2.1
* Non-maintainer upload.
* Add upstream patch fixing a crash in the open dialog
  in the zh_CN.utf8 locale. Closes: #487623.
  Thanks to Luca Bruno for the patch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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 :