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 Coord vtemp[order+1][order+1];
49 /* Copy control points */
50 std::copy(v, v+order+1, vtemp[0]);
52 /* Triangle computation */
53 for (unsigned i = 1; i <= order; i++) {
54 for (unsigned j = 0; j <= order - i; j++) {
55 vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]);
56 }
57 }
58 if(left != NULL)
59 for (unsigned j = 0; j <= order; j++)
60 left[j] = vtemp[j][0];
61 if(right != NULL)
62 for (unsigned j = 0; j <= order; j++)
63 right[j] = vtemp[order-j][j];
65 return (vtemp[order][0]);
66 }
69 class Bezier {
70 private:
71 std::valarray<Coord> c_;
73 friend Bezier portion(const Bezier & a, Coord from, Coord to);
75 friend OptInterval bounds_fast(Bezier const & b);
77 friend Bezier derivative(const Bezier & a);
79 protected:
80 Bezier(Coord const c[], unsigned ord) : c_(c, ord+1){
81 //std::copy(c, c+order()+1, &c_[0]);
82 }
84 public:
85 unsigned int order() const { return c_.size()-1;}
86 unsigned int size() const { return c_.size();}
88 Bezier() {}
89 Bezier(const Bezier& b) :c_(b.c_) {}
90 Bezier &operator=(Bezier const &other) {
91 if ( c_.size() != other.c_.size() ) {
92 c_.resize(other.c_.size());
93 }
94 c_ = other.c_;
95 return *this;
96 }
98 struct Order {
99 unsigned order;
100 explicit Order(Bezier const &b) : order(b.order()) {}
101 explicit Order(unsigned o) : order(o) {}
102 operator unsigned() const { return order; }
103 };
105 //Construct an arbitrary order bezier
106 Bezier(Order ord) : c_(0., ord.order+1) {
107 assert(ord.order == order());
108 }
110 explicit Bezier(Coord c0) : c_(0., 1) {
111 c_[0] = c0;
112 }
114 //Construct an order-1 bezier (linear Bézier)
115 Bezier(Coord c0, Coord c1) : c_(0., 2) {
116 c_[0] = c0; c_[1] = c1;
117 }
119 //Construct an order-2 bezier (quadratic Bézier)
120 Bezier(Coord c0, Coord c1, Coord c2) : c_(0., 3) {
121 c_[0] = c0; c_[1] = c1; c_[2] = c2;
122 }
124 //Construct an order-3 bezier (cubic Bézier)
125 Bezier(Coord c0, Coord c1, Coord c2, Coord c3) : c_(0., 4) {
126 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
127 }
129 void resize (unsigned int n, Coord v = 0)
130 {
131 c_.resize (n, v);
132 }
134 void clear()
135 {
136 c_.resize(0);
137 }
139 inline unsigned degree() const { return order(); }
141 //IMPL: FragmentConcept
142 typedef Coord output_type;
143 inline bool isZero() const {
144 for(unsigned i = 0; i <= order(); i++) {
145 if(c_[i] != 0) return false;
146 }
147 return true;
148 }
149 inline bool isConstant() const {
150 for(unsigned i = 1; i <= order(); i++) {
151 if(c_[i] != c_[0]) return false;
152 }
153 return true;
154 }
155 inline bool isFinite() const {
156 for(unsigned i = 0; i <= order(); i++) {
157 if(!IS_FINITE(c_[i])) return false;
158 }
159 return true;
160 }
161 inline Coord at0() const { return c_[0]; }
162 inline Coord at1() const { return c_[order()]; }
164 inline Coord valueAt(double t) const {
165 return subdivideArr(t, &c_[0], NULL, NULL, order());
166 }
167 inline Coord operator()(double t) const { return valueAt(t); }
169 SBasis toSBasis() const;
170 // inline SBasis toSBasis() const {
171 // SBasis sb;
172 // bezier_to_sbasis(sb, (*this));
173 // return sb;
174 // //return bezier_to_sbasis(&c_[0], order());
175 // }
177 //Only mutator
178 inline Coord &operator[](unsigned ix) { return c_[ix]; }
179 inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
180 inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
182 /* This is inelegant, as it uses several extra stores. I think there might be a way to
183 * evaluate roughly in situ. */
185 std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
186 std::vector<Coord> val_n_der;
187 Coord d_[order()+1];
188 unsigned nn = n_derivs + 1; // the size of the result vector equals n_derivs+1 ...
189 if(nn > order())
190 nn = order()+1; // .. but with a maximum of order() + 1!
191 for(unsigned i = 0; i < size(); i++)
192 d_[i] = c_[i];
193 for(unsigned di = 0; di < nn; di++) {
194 val_n_der.push_back(subdivideArr(t, d_, NULL, NULL, order() - di));
195 for(unsigned i = 0; i < order() - di; i++) {
196 d_[i] = (order()-di)*(d_[i+1] - d_[i]);
197 }
198 }
199 val_n_der.resize(n_derivs);
200 return val_n_der;
201 }
203 std::pair<Bezier, Bezier > subdivide(Coord t) const {
204 Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this));
205 subdivideArr(t, &c_[0], &a.c_[0], &b.c_[0], order());
206 return std::pair<Bezier, Bezier >(a, b);
207 }
209 std::vector<double> roots() const {
210 std::vector<double> solutions;
211 find_bernstein_roots(&c_[0], order(), solutions, 0, 0.0, 1.0);
212 return solutions;
213 }
214 };
217 void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
220 inline
221 SBasis Bezier::toSBasis() const {
222 SBasis sb;
223 bezier_to_sbasis(sb, (*this));
224 return sb;
225 //return bezier_to_sbasis(&c_[0], order());
226 }
228 //TODO: implement others
229 inline Bezier operator+(const Bezier & a, double v) {
230 Bezier result = Bezier(Bezier::Order(a));
231 for(unsigned i = 0; i <= a.order(); i++)
232 result[i] = a[i] + v;
233 return result;
234 }
236 inline Bezier operator-(const Bezier & a, double v) {
237 Bezier result = Bezier(Bezier::Order(a));
238 for(unsigned i = 0; i <= a.order(); i++)
239 result[i] = a[i] - v;
240 return result;
241 }
243 inline Bezier operator*(const Bezier & a, double v) {
244 Bezier result = Bezier(Bezier::Order(a));
245 for(unsigned i = 0; i <= a.order(); i++)
246 result[i] = a[i] * v;
247 return result;
248 }
250 inline Bezier operator/(const Bezier & a, double v) {
251 Bezier result = Bezier(Bezier::Order(a));
252 for(unsigned i = 0; i <= a.order(); i++)
253 result[i] = a[i] / v;
254 return result;
255 }
257 inline Bezier reverse(const Bezier & a) {
258 Bezier result = Bezier(Bezier::Order(a));
259 for(unsigned i = 0; i <= a.order(); i++)
260 result[i] = a[a.order() - i];
261 return result;
262 }
264 inline Bezier portion(const Bezier & a, double from, double to) {
265 //TODO: implement better?
266 Coord res[a.order()+1];
267 if(from == 0) {
268 if(to == 1) { return Bezier(a); }
269 subdivideArr(to, &a.c_[0], res, NULL, a.order());
270 return Bezier(res, a.order());
271 }
272 subdivideArr(from, &a.c_[0], NULL, res, a.order());
273 if(to == 1) return Bezier(res, a.order());
274 Coord res2[a.order()+1];
275 subdivideArr((to - from)/(1 - from), res, res2, NULL, a.order());
276 return Bezier(res2, a.order());
277 }
279 // XXX Todo: how to handle differing orders
280 inline std::vector<Point> bezier_points(const D2<Bezier > & a) {
281 std::vector<Point> result;
282 for(unsigned i = 0; i <= a[0].order(); i++) {
283 Point p;
284 for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
285 result.push_back(p);
286 }
287 return result;
288 }
290 inline Bezier derivative(const Bezier & a) {
291 //if(a.order() == 1) return Bezier(0.0);
292 if(a.order() == 1) return Bezier(a.c_[1]-a.c_[0]);
293 Bezier der(Bezier::Order(a.order()-1));
295 for(unsigned i = 0; i < a.order(); i++) {
296 der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]);
297 }
298 return der;
299 }
301 inline Bezier integral(const Bezier & a) {
302 Bezier inte(Bezier::Order(a.order()+1));
304 inte[0] = 0;
305 for(unsigned i = 0; i < inte.order(); i++) {
306 inte[i+1] = inte[i] + a[i]/(inte.order());
307 }
308 return inte;
309 }
311 inline OptInterval bounds_fast(Bezier const & b) {
312 return Interval::fromArray(&b.c_[0], b.size());
313 }
315 //TODO: better bounds exact
316 inline OptInterval bounds_exact(Bezier const & b) {
317 return bounds_exact(b.toSBasis());
318 }
320 inline OptInterval bounds_local(Bezier const & b, OptInterval i) {
321 //return bounds_local(b.toSBasis(), i);
322 if (i) {
323 return bounds_fast(portion(b, i->min(), i->max()));
324 } else {
325 return OptInterval();
326 }
327 }
329 inline std::ostream &operator<< (std::ostream &out_file, const Bezier & b) {
330 for(unsigned i = 0; i < b.size(); i++) {
331 out_file << b[i] << ", ";
332 }
333 return out_file;
334 }
336 }
337 #endif //SEEN_BEZIER_H
339 /*
340 Local Variables:
341 mode:c++
342 c-file-style:"stroustrup"
343 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
344 indent-tabs-mode:nil
345 fill-column:99
346 End:
347 */
348 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :