1fe846935eb21bb39acba881b9b14c134c76812a
1 /**
2 * \file bezier.h
3 * \brief \todo brief description
4 *
5 * Copyright 2007 MenTaLguY <mental@rydia.net>
6 * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
7 * Copyright 2007 Nathan Hurst <njh@njhurst.com>
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it either under the terms of the GNU Lesser General Public
11 * License version 2.1 as published by the Free Software Foundation
12 * (the "LGPL") or, at your option, under the terms of the Mozilla
13 * Public License Version 1.1 (the "MPL"). If you do not alter this
14 * notice, a recipient may use your version of this file under either
15 * the MPL or the LGPL.
16 *
17 * You should have received a copy of the LGPL along with this library
18 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 * You should have received a copy of the MPL along with this library
21 * in the file COPYING-MPL-1.1
22 *
23 * The contents of this file are subject to the Mozilla Public License
24 * Version 1.1 (the "License"); you may not use this file except in
25 * compliance with the License. You may obtain a copy of the License at
26 * http://www.mozilla.org/MPL/
27 *
28 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
30 * the specific language governing rights and limitations.
31 *
32 */
34 #ifndef SEEN_BEZIER_H
35 #define SEEN_BEZIER_H
37 #include <2geom/coord.h>
38 #include <valarray>
39 #include <2geom/isnan.h>
40 #include <2geom/d2.h>
41 #include <2geom/solver.h>
42 #include <boost/optional/optional.hpp>
44 namespace Geom {
46 inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, unsigned order) {
47 /*
48 * Bernstein :
49 * Evaluate a Bernstein function at a particular parameter value
50 * Fill in control points for resulting sub-curves.
51 *
52 */
54 unsigned N = order+1;
55 std::valarray<Coord> row(N);
56 for (unsigned i = 0; i < N; i++)
57 row[i] = v[i];
59 // Triangle computation
60 const double omt = (1-t);
61 if(left)
62 left[0] = row[0];
63 if(right)
64 right[order] = row[order];
65 for (unsigned i = 1; i < N; i++) {
66 for (unsigned j = 0; j < N - i; j++) {
67 row[j] = omt*row[j] + t*row[j+1];
68 }
69 if(left)
70 left[i] = row[0];
71 if(right)
72 right[order-i] = row[order-i];
73 }
74 return (row[0]);
75 /*
76 Coord vtemp[order+1][order+1];
78 // Copy control points
79 std::copy(v, v+order+1, vtemp[0]);
81 // Triangle computation
82 for (unsigned i = 1; i <= order; i++) {
83 for (unsigned j = 0; j <= order - i; j++) {
84 vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]);
85 }
86 }
87 if(left != NULL)
88 for (unsigned j = 0; j <= order; j++)
89 left[j] = vtemp[j][0];
90 if(right != NULL)
91 for (unsigned j = 0; j <= order; j++)
92 right[j] = vtemp[order-j][j];
94 return (vtemp[order][0]);*/
95 }
97 template <typename T>
98 inline T bernsteinValueAt(double t, T const *c_, unsigned n) {
99 double u = 1.0 - t;
100 double bc = 1;
101 double tn = 1;
102 T tmp = c_[0]*u;
103 for(unsigned i=1; i<n; i++){
104 tn = tn*t;
105 bc = bc*(n-i+1)/i;
106 tmp = (tmp + tn*bc*c_[i])*u;
107 }
108 return (tmp + tn*t*c_[n]);
109 }
112 class Bezier {
113 private:
114 std::valarray<Coord> c_;
116 friend Bezier portion(const Bezier & a, Coord from, Coord to);
118 friend OptInterval bounds_fast(Bezier const & b);
120 friend Bezier derivative(const Bezier & a);
122 protected:
123 Bezier(Coord const c[], unsigned ord) : c_(c, ord+1){
124 //std::copy(c, c+order()+1, &c_[0]);
125 }
127 public:
128 unsigned int order() const { return c_.size()-1;}
129 unsigned int size() const { return c_.size();}
131 Bezier() {}
132 Bezier(const Bezier& b) :c_(b.c_) {}
133 Bezier &operator=(Bezier const &other) {
134 if ( c_.size() != other.c_.size() ) {
135 c_.resize(other.c_.size());
136 }
137 c_ = other.c_;
138 return *this;
139 }
141 struct Order {
142 unsigned order;
143 explicit Order(Bezier const &b) : order(b.order()) {}
144 explicit Order(unsigned o) : order(o) {}
145 operator unsigned() const { return order; }
146 };
148 //Construct an arbitrary order bezier
149 Bezier(Order ord) : c_(0., ord.order+1) {
150 assert(ord.order == order());
151 }
153 explicit Bezier(Coord c0) : c_(0., 1) {
154 c_[0] = c0;
155 }
157 //Construct an order-1 bezier (linear Bézier)
158 Bezier(Coord c0, Coord c1) : c_(0., 2) {
159 c_[0] = c0; c_[1] = c1;
160 }
162 //Construct an order-2 bezier (quadratic Bézier)
163 Bezier(Coord c0, Coord c1, Coord c2) : c_(0., 3) {
164 c_[0] = c0; c_[1] = c1; c_[2] = c2;
165 }
167 //Construct an order-3 bezier (cubic Bézier)
168 Bezier(Coord c0, Coord c1, Coord c2, Coord c3) : c_(0., 4) {
169 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
170 }
172 void resize (unsigned int n, Coord v = 0)
173 {
174 c_.resize (n, v);
175 }
177 void clear()
178 {
179 c_.resize(0);
180 }
182 inline unsigned degree() const { return order(); }
184 //IMPL: FragmentConcept
185 typedef Coord output_type;
186 inline bool isZero() const {
187 for(unsigned i = 0; i <= order(); i++) {
188 if(c_[i] != 0) return false;
189 }
190 return true;
191 }
192 inline bool isConstant() const {
193 for(unsigned i = 1; i <= order(); i++) {
194 if(c_[i] != c_[0]) return false;
195 }
196 return true;
197 }
198 inline bool isFinite() const {
199 for(unsigned i = 0; i <= order(); i++) {
200 if(!IS_FINITE(c_[i])) return false;
201 }
202 return true;
203 }
204 inline Coord at0() const { return c_[0]; }
205 inline Coord at1() const { return c_[order()]; }
207 inline Coord valueAt(double t) const {
208 int n = order();
209 double u, bc, tn, tmp;
210 int i;
211 u = 1.0 - t;
212 bc = 1;
213 tn = 1;
214 tmp = c_[0]*u;
215 for(i=1; i<n; i++){
216 tn = tn*t;
217 bc = bc*(n-i+1)/i;
218 tmp = (tmp + tn*bc*c_[i])*u;
219 }
220 return (tmp + tn*t*c_[n]);
221 //return subdivideArr(t, &c_[0], NULL, NULL, order());
222 }
223 inline Coord operator()(double t) const { return valueAt(t); }
225 SBasis toSBasis() const;
226 // inline SBasis toSBasis() const {
227 // SBasis sb;
228 // bezier_to_sbasis(sb, (*this));
229 // return sb;
230 // //return bezier_to_sbasis(&c_[0], order());
231 // }
233 //Only mutator
234 inline Coord &operator[](unsigned ix) { return c_[ix]; }
235 inline Coord const &operator[](unsigned ix) const { return const_cast<std::valarray<Coord>&>(c_)[ix]; }
236 //inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
237 inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
239 /* This is inelegant, as it uses several extra stores. I think there might be a way to
240 * evaluate roughly in situ. */
242 std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
243 std::vector<Coord> val_n_der;
244 std::valarray<Coord> d_(order()+1);
245 unsigned nn = n_derivs + 1; // the size of the result vector equals n_derivs+1 ...
246 if(nn > order())
247 nn = order()+1; // .. but with a maximum of order() + 1!
248 for(unsigned i = 0; i < size(); i++)
249 d_[i] = c_[i];
250 val_n_der.resize(nn);
251 for(unsigned di = 0; di < nn; di++) {
252 //val_n_der[di] = (subdivideArr(t, &d_[0], NULL, NULL, order() - di));
253 val_n_der[di] = bernsteinValueAt(t, &d_[0], order() - di);
254 for(unsigned i = 0; i < order() - di; i++) {
255 d_[i] = (order()-di)*(d_[i+1] - d_[i]);
256 }
257 }
258 return val_n_der;
259 }
261 std::pair<Bezier, Bezier > subdivide(Coord t) const {
262 Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this));
263 subdivideArr(t, &const_cast<std::valarray<Coord>&>(c_)[0], &a.c_[0], &b.c_[0], order());
264 return std::pair<Bezier, Bezier >(a, b);
265 }
267 std::vector<double> roots() const {
268 std::vector<double> solutions;
269 find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, 0.0, 1.0);
270 return solutions;
271 }
272 std::vector<double> roots(Interval const ivl) const {
273 std::vector<double> solutions;
274 find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl[0], ivl[1]);
275 return solutions;
276 }
277 };
280 void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
283 inline
284 SBasis Bezier::toSBasis() const {
285 SBasis sb;
286 bezier_to_sbasis(sb, (*this));
287 return sb;
288 //return bezier_to_sbasis(&c_[0], order());
289 }
291 //TODO: implement others
292 inline Bezier operator+(const Bezier & a, double v) {
293 Bezier result = Bezier(Bezier::Order(a));
294 for(unsigned i = 0; i <= a.order(); i++)
295 result[i] = a[i] + v;
296 return result;
297 }
299 inline Bezier operator-(const Bezier & a, double v) {
300 Bezier result = Bezier(Bezier::Order(a));
301 for(unsigned i = 0; i <= a.order(); i++)
302 result[i] = a[i] - v;
303 return result;
304 }
306 inline Bezier operator*(const Bezier & a, double v) {
307 Bezier result = Bezier(Bezier::Order(a));
308 for(unsigned i = 0; i <= a.order(); i++)
309 result[i] = a[i] * v;
310 return result;
311 }
313 inline Bezier operator/(const Bezier & a, double v) {
314 Bezier result = Bezier(Bezier::Order(a));
315 for(unsigned i = 0; i <= a.order(); i++)
316 result[i] = a[i] / v;
317 return result;
318 }
320 inline Bezier reverse(const Bezier & a) {
321 Bezier result = Bezier(Bezier::Order(a));
322 for(unsigned i = 0; i <= a.order(); i++)
323 result[i] = a[a.order() - i];
324 return result;
325 }
327 inline Bezier portion(const Bezier & a, double from, double to) {
328 //TODO: implement better?
329 std::valarray<Coord> res(a.order() + 1);
330 if(from == 0) {
331 if(to == 1) { return Bezier(a); }
332 subdivideArr(to, &const_cast<Bezier&>(a).c_[0], &res[0], NULL, a.order());
333 return Bezier(&res[0], a.order());
334 }
335 subdivideArr(from, &const_cast<Bezier&>(a).c_[0], NULL, &res[0], a.order());
336 if(to == 1) return Bezier(&res[0], a.order());
337 std::valarray<Coord> res2(a.order()+1);
338 subdivideArr((to - from)/(1 - from), &res[0], &res2[0], NULL, a.order());
339 return Bezier(&res2[0], a.order());
340 }
342 // XXX Todo: how to handle differing orders
343 inline std::vector<Point> bezier_points(const D2<Bezier > & a) {
344 std::vector<Point> result;
345 for(unsigned i = 0; i <= a[0].order(); i++) {
346 Point p;
347 for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
348 result.push_back(p);
349 }
350 return result;
351 }
353 inline Bezier derivative(const Bezier & a) {
354 //if(a.order() == 1) return Bezier(0.0);
355 if(a.order() == 1) return Bezier(a.c_[1]-a.c_[0]);
356 Bezier der(Bezier::Order(a.order()-1));
358 for(unsigned i = 0; i < a.order(); i++) {
359 der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]);
360 }
361 return der;
362 }
364 inline Bezier integral(const Bezier & a) {
365 Bezier inte(Bezier::Order(a.order()+1));
367 inte[0] = 0;
368 for(unsigned i = 0; i < inte.order(); i++) {
369 inte[i+1] = inte[i] + a[i]/(inte.order());
370 }
371 return inte;
372 }
374 inline OptInterval bounds_fast(Bezier const & b) {
375 return Interval::fromArray(&const_cast<Bezier&>(b).c_[0], b.size());
376 }
378 //TODO: better bounds exact
379 inline OptInterval bounds_exact(Bezier const & b) {
380 return bounds_exact(b.toSBasis());
381 }
383 inline OptInterval bounds_local(Bezier const & b, OptInterval i) {
384 //return bounds_local(b.toSBasis(), i);
385 if (i) {
386 return bounds_fast(portion(b, i->min(), i->max()));
387 } else {
388 return OptInterval();
389 }
390 }
392 inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
393 for(unsigned i = 0; i < b.size(); i++) {
394 out_file << b[i] << ", ";
395 }
396 return out_file;
397 }
399 }
400 #endif //SEEN_BEZIER_H
402 /*
403 Local Variables:
404 mode:c++
405 c-file-style:"stroustrup"
406 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
407 indent-tabs-mode:nil
408 fill-column:99
409 End:
410 */
411 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :