Inkscape
Vector Graphics Editor
Loading...
Searching...
No Matches
sbasis-of.h
Go to the documentation of this file.
1
36#ifndef SEEN_SBASIS_OF_H
37#define SEEN_SBASIS_OF_H
38#include <vector>
39#include <cassert>
40#include <iostream>
41
42#include <2geom/interval.h>
43#include <2geom/utils.h>
44#include <2geom/exception.h>
45
47
48namespace Geom {
49
50template<typename T>
51class SBasisOf;
52
53#ifdef USE_SBASIS_OF
55#endif
56
57/*** An empty SBasisOf<T> is identically 0. */
58template<typename T>
59class SBasisOf : public std::vector<LinearOf<T> >{
60public:
62 explicit SBasisOf(T a) {
63 this->push_back(LinearOf<T>(a,a));
64 }
65 SBasisOf(SBasisOf<T> const & a) :
66 std::vector<LinearOf<T> >(a)
67 {}
68 SBasisOf(LinearOf<T> const & bo) {
69 this->push_back(bo);
70 }
72 this->push_back(*bo);
73 }
74 //static unsigned input_dim(){return T::input_dim()+1;}
75
76 //IMPL: FragmentConcept
77 typedef T output_type;
78 inline bool isZero() const {
79 if(this->empty()) return true;
80 for(unsigned i = 0; i < this->size(); i++) {
81 if(!(*this)[i].isZero()) return false;
82 }
83 return true;
84 }
85 inline bool isConstant() const {
86 if (this->empty()) return true;
87 for (unsigned i = 0; i < this->size(); i++) {
88 if(!(*this)[i].isConstant()) return false;
89 }
90 return true;
91 }
92
93 //TODO: code this...
94 bool isFinite() const;
95
96 inline T at0() const {
97 if(this->empty()) return T(0); else return (*this)[0][0];
98 }
99 inline T at1() const{
100 if(this->empty()) return T(0); else return (*this)[0][1];
101 }
102
103 T valueAt(double t) const {
104 double s = t*(1-t);
105 T p0 = T(0.), p1 = T(0.);
106 for(unsigned k = this->size(); k > 0; k--) {
107 const LinearOf<T> &lin = (*this)[k-1];
108 p0 = p0*s + lin[0];
109 p1 = p1*s + lin[1];
110 }
111 return p0*(1-t) + p1*t;
112 }
113
114 T operator()(double t) const {
115 return valueAt(t);
116 }
117
121 std::vector<T> valueAndDerivatives(double t, unsigned n) const{
122 std::vector<T> ret(n+1);
123 ret[0] = valueAt(t);
124 SBasisOf<T> tmp = *this;
125 for(unsigned i = 0; i < n; i++) {
126 tmp.derive();
127 ret[i+1] = tmp.valueAt(t);
128 }
129 return ret;
130 }
131
132 //The following lines only makes sens if T=double!
133 SBasisOf<T> toSBasis() const { return SBasisOf<T>(*this); }
134 double tailError(unsigned tail) const{
135 Interval bs = *bounds_fast(*this, tail);
136 return std::max(fabs(bs.min()),fabs(bs.max()));
137 }
138
139// compute f(g)
141
142 LinearOf<T> operator[](unsigned i) const {
143 assert(i < this->size());
144 return std::vector<LinearOf<T> >::operator[](i);
145 }
146
147//MUTATOR PRISON
148 LinearOf<T>& operator[](unsigned i) { return this->at(i); }
149
150 //remove extra zeros
151 void normalize() {
152 while(!this->empty() && this->back().isZero())
153 this->pop_back();
154 }
155
156 void truncate(unsigned k) { if(k < this->size()) this->resize(k); }
157private:
158 void derive(); // in place version
159 unsigned dim;
160};
161
162//template<>
163//inline unsigned SBasisOf<double>::input_dim() { return 1; }
164
165//--------------------------------------------------------------------------
166#ifdef USE_SBASIS_OF
167
168//implemented in sbasis-roots.cpp
171OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0);
172
173std::vector<double> roots(SBasis const & s);
174std::vector<std::vector<double> > multi_roots(SBasis const &f,
175 std::vector<double> const &levels,
176 double htol=1e-7,
177 double vtol=1e-7,
178 double a=0,
179 double b=1);
180#endif
181//--------------------------------------------------------------------------
182
183
184//TODO: figure out how to stick this in linear, while not adding an sbasis dep
185template<typename T>
186inline SBasisOf<T> LinearOf<T>::toSBasis() const { return SBasisOf<T>(*this); }
187
188template<typename T>
191 result.reserve(a.size());
192 for(unsigned k = 0; k < a.size(); k++)
193 result.push_back(reverse(a[k]));
194 return result;
195}
196
197//IMPL: ScalableConcept
198template<typename T>
200 if(p.isZero()) return SBasisOf<T>();
202 result.reserve(p.size());
203
204 for(unsigned i = 0; i < p.size(); i++) {
205 result.push_back(-p[i]);
206 }
207 return result;
208}
209
210template<typename T>
213 //TODO: what does this mean for vectors of vectors??
214 //c.reserve(a.size());
215 for(unsigned i = 0; i < a.size(); i++)
216 c.push_back(a[i] * k);
217 return c;
218}
219
220template<typename T>
221inline SBasisOf<T> operator*(double k, SBasisOf<T> const &a) { return a*k; }
222template<typename T>
223inline SBasisOf<T> operator/(SBasisOf<T> const &a, double k) { return a*(1./k); }
224template<typename T>
226 if (a.isZero()) return a;
227 if (b == 0)
228 a.clear();
229 else
230 for(unsigned i = 0; i < a.size(); i++)
231 a[i] *= b;
232 return a;
233}
234
235template<typename T>
236inline SBasisOf<T>& operator/=(SBasisOf<T>& a, double b) { return (a*=(1./b)); }
237
238/*
239//We can also multiply by element of ring coeff T:
240template<typename T>
241SBasisOf<T> operator*(SBasisOf<T> const &a, T k){
242 SBasisOf<T> c;
243 //TODO: what does this mean for vectors of vectors??
244 //c.reserve(a.size());
245 for(unsigned i = 0; i < a.size(); i++)
246 c.push_back(a[i] * k);
247 return c;
248}
249
250template<typename T>
251inline SBasisOf<T> operator*(T k, SBasisOf<T> const &a) { return a*k; }
252template<typename T>
253SBasisOf<T>& operator*=(SBasisOf<T>& a, T b){
254 if (a.isZero()) return a;
255 if (b == 0)
256 a.clear();
257 else
258 for(unsigned i = 0; i < a.size(); i++)
259 a[i] *= b;
260 return a;
261}
262*/
263
264//IMPL: AddableConcept
265template<typename T>
266inline SBasisOf<T> operator+(const SBasisOf<T>& a, const SBasisOf<T>& b){
268 const unsigned out_size = std::max(a.size(), b.size());
269 const unsigned min_size = std::min(a.size(), b.size());
270 //TODO: what does this mean for vector<vector>;
271 //result.reserve(out_size);
272
273 for(unsigned i = 0; i < min_size; i++) {
274 result.push_back(a[i] + b[i]);
275 }
276 for(unsigned i = min_size; i < a.size(); i++)
277 result.push_back(a[i]);
278 for(unsigned i = min_size; i < b.size(); i++)
279 result.push_back(b[i]);
280
281 assert(result.size() == out_size);
282 return result;
283}
284
285template<typename T>
288 const unsigned out_size = std::max(a.size(), b.size());
289 const unsigned min_size = std::min(a.size(), b.size());
290 //TODO: what does this mean for vector<vector>;
291 //result.reserve(out_size);
292
293 for(unsigned i = 0; i < min_size; i++) {
294 result.push_back(a[i] - b[i]);
295 }
296 for(unsigned i = min_size; i < a.size(); i++)
297 result.push_back(a[i]);
298 for(unsigned i = min_size; i < b.size(); i++)
299 result.push_back(-b[i]);
300
301 assert(result.size() == out_size);
302 return result;
303}
304
305template<typename T>
307 const unsigned out_size = std::max(a.size(), b.size());
308 const unsigned min_size = std::min(a.size(), b.size());
309 //TODO: what does this mean for vectors of vectors
310 //a.reserve(out_size);
311 for(unsigned i = 0; i < min_size; i++)
312 a[i] += b[i];
313 for(unsigned i = min_size; i < b.size(); i++)
314 a.push_back(b[i]);
315
316 assert(a.size() == out_size);
317 return a;
318}
319
320template<typename T>
322 const unsigned out_size = std::max(a.size(), b.size());
323 const unsigned min_size = std::min(a.size(), b.size());
324 //TODO: what does this mean for vectors of vectors
325 //a.reserve(out_size);
326 for(unsigned i = 0; i < min_size; i++)
327 a[i] -= b[i];
328 for(unsigned i = min_size; i < b.size(); i++)
329 a.push_back(-b[i]);
330
331 assert(a.size() == out_size);
332 return a;
333}
334
335//TODO: remove?
336template<typename T>
337inline SBasisOf<T> operator+(const SBasisOf<T> & a, LinearOf<T> const & b) {
338 if(b.isZero()) return a;
339 if(a.isZero()) return b;
341 result[0] += b;
342 return result;
343}
344template<typename T>
345inline SBasisOf<T> operator-(const SBasisOf<T> & a, LinearOf<T> const & b) {
346 if(b.isZero()) return a;
348 result[0] -= b;
349 return result;
350}
351template<typename T>
353 if(a.isZero())
354 a.push_back(b);
355 else
356 a[0] += b;
357 return a;
358}
359template<typename T>
361 if(a.isZero())
362 a.push_back(-b);
363 else
364 a[0] -= b;
365 return a;
366}
367
368//IMPL: OffsetableConcept
369/*
370template<typename T>
371inline SBasisOf<T> operator+(const SBasisOf<T> & a, double b) {
372 if(a.isZero()) return LinearOf<T>(b, b);
373 SBasisOf<T> result(a);
374 result[0] += b;
375 return result;
376}
377template<typename T>
378inline SBasisOf<T> operator-(const SBasisOf<T> & a, double b) {
379 if(a.isZero()) return LinearOf<T>(-b, -b);
380 SBasisOf<T> result(a);
381 result[0] -= b;
382 return result;
383}
384template<typename T>
385inline SBasisOf<T>& operator+=(SBasisOf<T>& a, double b) {
386 if(a.isZero())
387 a.push_back(LinearOf<T>(b,b));
388 else
389 a[0] += b;
390 return a;
391}
392template<typename T>
393inline SBasisOf<T>& operator-=(SBasisOf<T>& a, double b) {
394 if(a.isZero())
395 a.push_back(LinearOf<T>(-b,-b));
396 else
397 a[0] -= b;
398 return a;
399}
400*/
401//We can also offset by elements of coeff ring T
402template<typename T>
403inline SBasisOf<T> operator+(const SBasisOf<T> & a, T b) {
404 if(a.isZero()) return LinearOf<T>(b, b);
406 result[0] += b;
407 return result;
408}
409template<typename T>
410inline SBasisOf<T> operator-(const SBasisOf<T> & a, T b) {
411 if(a.isZero()) return LinearOf<T>(-b, -b);
413 result[0] -= b;
414 return result;
415}
416template<typename T>
418 if(a.isZero())
419 a.push_back(LinearOf<T>(b,b));
420 else
421 a[0] += b;
422 return a;
423}
424template<typename T>
426 if(a.isZero())
427 a.push_back(LinearOf<T>(-b,-b));
428 else
429 a[0] -= b;
430 return a;
431}
432
433
434template<typename T>
435SBasisOf<T> shift(SBasisOf<T> const &a, int sh){
436 SBasisOf<T> c = a;
437 if(sh > 0) {
438 c.insert(c.begin(), sh, LinearOf<T>(0,0));
439 } else {
440 //TODO: truncate
441 }
442 return c;
443}
444
445template<typename T>
446SBasisOf<T> shift(LinearOf<T> const &a, int sh) {
448 if(sh > 0) {
449 c.insert(c.begin(), sh, LinearOf<T>(0,0));
450 c.push_back(a);
451 }
452 return c;
453}
454
455template<typename T>
456inline SBasisOf<T> truncate(SBasisOf<T> const &a, unsigned terms) {
458 c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));
459 return c;
460}
461
462template<typename T>
464 if(a.isZero() || b.isZero())
465 return c;
466 c.resize(a.size() + b.size(), LinearOf<T>(T(0.),T(0.)));
467 for(unsigned j = 0; j < b.size(); j++) {
468 for(unsigned i = j; i < a.size()+j; i++) {
469 T tri = (b[j][1]-b[j][0])*(a[i-j][1]-a[i-j][0]);
470 c[i+1/*shift*/] += LinearOf<T>(-tri);
471 }
472 }
473 for(unsigned j = 0; j < b.size(); j++) {
474 for(unsigned i = j; i < a.size()+j; i++) {
475 for(unsigned dim = 0; dim < 2; dim++)
476 c[i][dim] += b[j][dim]*a[i-j][dim];
477 }
478 }
479 c.normalize();
480 //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
481 return c;
482}
483
484template<typename T>
487 if(a.isZero() || b.isZero())
488 return c;
489 return multiply_add(a, b, c);
490}
491
492template<typename T>
494 SBasisOf<T> a;
495 T aTri = T(0.);
496 for(int k = c.size()-1; k >= 0; k--) {
497 aTri = (HatOf<T>(c[k]).d + (k+1)*aTri/2)/(2*k+1);
498 a[k][0] -= aTri/2;
499 a[k][1] += aTri/2;
500 }
501 a.normalize();
502 return a;
503}
504
505template<typename T>
508 c.resize(a.size(), LinearOf<T>());
509 if(a.isZero())
510 return c;
511
512 for(unsigned k = 0; k < a.size()-1; k++) {
513 T d = (2*k+1)*(a[k][1] - a[k][0]);
514
515 c[k][0] = d + (k+1)*a[k+1][0];
516 c[k][1] = d - (k+1)*a[k+1][1];
517 }
518 int k = a.size()-1;
519 T d = (2*k+1)*(a[k][1] - a[k][0]);
520 //TODO: do a real test to know if d==0!
521 if(d == T(0.0))
522 c.pop_back();
523 else {
524 c[k][0] = d;
525 c[k][1] = d;
526 }
527
528 return c;
529}
530
531template<typename T>
532void SBasisOf<T>::derive() { // in place version
533 if(isZero()) return;
534 for(unsigned k = 0; k < this->size()-1; k++) {
535 T d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
536
537 (*this)[k][0] = d + (k+1)*(*this)[k+1][0];
538 (*this)[k][1] = d - (k+1)*(*this)[k+1][1];
539 }
540 int k = this->size()-1;
541 T d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
542 if(d == 0)//TODO: give this a meaning for general coeff ring.
543 this->pop_back();
544 else {
545 (*this)[k][0] = d;
546 (*this)[k][1] = d;
547 }
548}
549
550
551template<typename T>
552inline SBasisOf<T> operator*(SBasisOf<T> const & a, SBasisOf<T> const & b) {
553 return multiply(a, b);
554}
555
556template<typename T>
558 a = multiply(a, b);
559 return a;
560}
561
562// a(b(t))
563//TODO: compose all compatibles types!
564template<typename T>
567 SBasisOf<T> r;
568
569 for(int i = a.size()-1; i >= 0; i--) {
570 r = multiply_add(r, s, SBasisOf<T>(LinearOf<T>(HatOf<T>(a[i][0]))) - b*a[i][0] + b*a[i][1]);
571 }
572 return r;
573}
574
575template<typename T>
576SBasisOf<T> compose(SBasisOf<T> const &a, SBasisOf<T> const &b, unsigned k){
578 SBasisOf<T> r;
579
580 for(int i = a.size()-1; i >= 0; i--) {
581 r = multiply_add(r, s, SBasisOf<T>(LinearOf<T>(HatOf<T>(a[i][0]))) - b*a[i][0] + b*a[i][1]);
582 }
583 r.truncate(k);
584 return r;
585}
586template<typename T>
588 return compose(SBasisOf<T>(a),b);
589}
590template<typename T>
592 return compose(a,SBasisOf<T>(b));
593}
594template<typename T>//TODO: don't be so lazy!!
596 return compose(SBasisOf<T>(a),SBasisOf<T>(b));
597}
598
599
600
601template<typename T>
602inline SBasisOf<T> portion(const SBasisOf<T> &t, double from, double to) { return compose(t, LinearOf<T>(from, to)); }
603
604// compute f(g)
605template<typename T>
606inline SBasisOf<T>
608 return compose(*this, g);
609}
610
611template<typename T>
612inline std::ostream &operator<< (std::ostream &out_file, const LinearOf<T> &bo) {
613 out_file << "{" << bo[0] << ", " << bo[1] << "}";
614 return out_file;
615}
616
617template<typename T>
618inline std::ostream &operator<< (std::ostream &out_file, const SBasisOf<T> & p) {
619 for(unsigned i = 0; i < p.size(); i++) {
620 out_file << p[i] << "s^" << i << " + ";
621 }
622 return out_file;
623}
624
625};
626#endif
627
628
629/*
630 Local Variables:
631 mode:c++
632 c-file-style:"stroustrup"
633 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
634 indent-tabs-mode:nil
635 fill-column:99
636 End:
637*/
638// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
Defines the different types of exceptions that 2geom can throw.
Various utility functions.
Range of real numbers that is never empty.
Definition interval.h:59
SBasisOf< T > toSBasis() const
Definition sbasis-of.h:186
bool isZero() const
Definition linear-of.h:107
Range of real numbers that can be empty.
Definition interval.h:199
SBasisOf< T > toSBasis() const
Definition sbasis-of.h:133
bool isConstant() const
Definition sbasis-of.h:85
T valueAt(double t) const
Definition sbasis-of.h:103
SBasisOf< T > operator()(SBasisOf< T > const &g) const
Definition sbasis-of.h:607
SBasisOf(SBasisOf< T > const &a)
Definition sbasis-of.h:65
LinearOf< T > operator[](unsigned i) const
Definition sbasis-of.h:142
T at1() const
Definition sbasis-of.h:99
bool isFinite() const
LinearOf< T > & operator[](unsigned i)
Definition sbasis-of.h:148
double tailError(unsigned tail) const
Definition sbasis-of.h:134
SBasisOf(LinearOf< T > *bo)
Definition sbasis-of.h:71
T operator()(double t) const
Definition sbasis-of.h:114
T at0() const
Definition sbasis-of.h:96
void truncate(unsigned k)
Definition sbasis-of.h:156
SBasisOf(LinearOf< T > const &bo)
Definition sbasis-of.h:68
unsigned dim
Definition sbasis-of.h:159
bool isZero() const
Definition sbasis-of.h:78
void normalize()
Definition sbasis-of.h:151
std::vector< T > valueAndDerivatives(double t, unsigned n) const
The size of the returned vector equals n+1.
Definition sbasis-of.h:121
Polynomial in symmetric power basis.
Definition sbasis.h:70
Css & result
Geom::IntPoint size
double c[8][4]
const unsigned order
Simple closed interval class.
Linear fragment function class.
Various utility functions.
Definition affine.h:22
OptInterval bounds_exact(Bezier const &b)
Definition bezier.cpp:310
Bezier reverse(const Bezier &a)
Definition bezier.h:342
OptInterval bounds_local(Bezier const &b, OptInterval const &i)
Definition bezier.cpp:320
Bezier multiply(Bezier const &f, Bezier const &g)
Definition bezier.h:337
D2< T > operator+=(D2< T > &a, D2< T > const &b)
Definition d2.h:219
D2< T > operator/(D2< T > const &a, Point const &b)
Definition d2.h:258
D2< T > operator-(D2< T > const &a, D2< T > const &b)
Definition d2.h:209
D2< T > operator*=(D2< T > &a, Point const &b)
Definition d2.h:268
SBasisOf< T > shift(SBasisOf< T > const &a, int sh)
Definition sbasis-of.h:435
Bezier operator*(Bezier const &f, Bezier const &g)
Definition bezier.cpp:221
D2< T > operator-=(D2< T > &a, D2< T > const &b)
Definition d2.h:228
D2< T > compose(D2< T > const &a, T const &b)
Definition d2.h:405
std::vector< double > roots(SBasis const &s)
Bezier portion(const Bezier &a, double from, double to)
Definition bezier.cpp:250
SBasisOf< T > multiply_add(SBasisOf< T > const &a, SBasisOf< T > const &b, SBasisOf< T > c)
Definition sbasis-of.h:463
Bezier integral(Bezier const &a)
Definition bezier.cpp:294
Bezier derivative(Bezier const &a)
Definition bezier.cpp:282
D2< T > operator+(D2< T > const &a, D2< T > const &b)
Definition d2.h:199
std::vector< std::vector< double > > multi_roots(SBasis const &f, std::vector< double > const &levels, double htol=1e-7, double vtol=1e-7, double a=0, double b=1)
D2< SBasis > truncate(D2< SBasis > const &a, unsigned terms)
Definition d2-sbasis.cpp:52
OptInterval bounds_fast(Bezier const &b)
Definition bezier.cpp:305
D2< T > operator/=(D2< T > &a, Point const &b)
Definition d2.h:277
STL namespace.
unsigned long bs
Definition quantize.cpp:38
void resize()
Definition resize.cpp:80