31 #ifndef SEEN_GEOM_PW_SB_H
32 #define SEEN_GEOM_PW_SB_H
40 #include <boost/concept_check.hpp>
49 std::vector<double> cuts;
61 typedef typename T::output_type output_type;
63 explicit Piecewise(
const output_type & v) {
69 inline T operator[](
unsigned i)
const {
return segs[i]; }
70 inline T &operator[](
unsigned i) {
return segs[i]; }
71 inline output_type operator()(
double t)
const {
return valueAt(t); }
72 inline output_type valueAt(
double t)
const {
74 return segs[n](
segT(t, n));
80 inline unsigned size()
const {
return segs.size(); }
81 inline bool empty()
const {
return segs.empty(); }
86 inline void push(
const T &s,
double to) {
87 assert(cuts.size() - segs.size() == 1);
92 inline void push_cut(
double c) {
93 assert_invariants(cuts.empty() || c > cuts.back());
97 inline void push_seg(
const T &s) { segs.push_back(s); }
102 inline unsigned segN(
double t,
int low = 0,
int high = -1)
const {
103 high = (high == -1) ? size() : high;
104 if(t < cuts[0])
return 0;
105 if(t >= cuts[size()])
return size() - 1;
107 int mid = (high + low) / 2;
108 double mv = cuts[mid];
110 if(t < cuts[mid + 1])
return mid;
else low = mid + 1;
112 if(cuts[mid - 1] < t)
return mid - 1;
else high = mid - 1;
124 inline double segT(
double t,
int i = -1)
const {
125 if(i == -1) i =
segN(t);
127 return (t - cuts[i]) / (cuts[i+1] - cuts[i]);
130 inline double mapToDomain(
double t,
unsigned i)
const {
131 return (1-t)*cuts[i] + t*cuts[i+1];
135 inline void offsetDomain(
double o) {
136 assert(is_finite(o));
138 for(
unsigned i = 0; i <= size(); i++)
143 inline void scaleDomain(
double s) {
146 cuts.clear(); segs.clear();
149 for(
unsigned i = 0; i <= size(); i++)
154 inline Interval domain()
const {
return Interval(cuts.front(), cuts.back()); }
157 inline void setDomain(Interval dom) {
160 cuts.clear(); segs.clear();
163 double cf = cuts.front();
164 double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
165 for(
unsigned i = 0; i <= size(); i++)
166 cuts[i] = (cuts[i] - cf) * s + o;
170 inline void concat(
const Piecewise<T> &other) {
171 if(other.empty())
return;
174 cuts = other.cuts; segs = other.segs;
178 segs.insert(segs.end(), other.segs.begin(), other.segs.end());
179 double t = cuts.back() - other.cuts.front();
180 for(
unsigned i = 0; i < other.size(); i++)
181 push_cut(other.cuts[i + 1] + t);
185 inline void continuousConcat(
const Piecewise<T> &other) {
186 boost::function_requires<AddableConcept<typename T::output_type> >();
187 if(other.empty())
return;
188 typename T::output_type y = segs.back().at1() - other.segs.front().at0();
191 for(
unsigned i = 0; i < other.size(); i++)
192 push_seg(other[i] + y);
197 double t = cuts.back() - other.cuts.front();
198 for(
unsigned i = 0; i < other.size(); i++)
199 push(other[i] + y, other.cuts[i + 1] + t);
203 inline bool invariants()
const {
205 if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))
208 for(
unsigned i = 0; i < segs.size(); i++)
209 if(cuts[i] >= cuts[i+1])
217 inline typename FragmentConcept<T>::BoundsType bounds_fast(
const Piecewise<T> &f) {
218 boost::function_requires<FragmentConcept<T> >();
220 if(f.empty())
return typename FragmentConcept<T>::BoundsType();
221 typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));
222 for(
unsigned i = 1; i < f.size(); i++)
223 ret.unionWith(bounds_fast(f[i]));
228 inline typename FragmentConcept<T>::BoundsType bounds_exact(
const Piecewise<T> &f) {
229 boost::function_requires<FragmentConcept<T> >();
231 if(f.empty())
return typename FragmentConcept<T>::BoundsType();
232 typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));
233 for(
unsigned i = 1; i < f.size(); i++)
234 ret.unionWith(bounds_exact(f[i]));
239 inline typename FragmentConcept<T>::BoundsType bounds_local(
const Piecewise<T> &f,
const Interval &m) {
240 boost::function_requires<FragmentConcept<T> >();
242 if(f.empty())
return typename FragmentConcept<T>::BoundsType();
243 if(m.isEmpty())
return typename FragmentConcept<T>::BoundsType(f(m.min()));
245 unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
246 double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
248 if(fi == ti)
return bounds_local(f[fi], Interval(ft, tt));
250 typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));
251 for(
unsigned i = fi + 1; i < ti; i++)
252 ret.unionWith(bounds_exact(f[i]));
253 if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));
260 T elem_portion(
const Piecewise<T> &a,
unsigned i,
double from,
double to) {
261 assert(i < a.size());
262 double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);
263 return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );
277 assert(pw.invariants());
281 ret.cuts.reserve(c.size() + pw.cuts.size());
282 ret.segs.reserve(c.size() + pw.cuts.size() - 1);
286 for(
unsigned i = 0; i < c.size() - 1; i++)
291 unsigned si = 0, ci = 0;
294 while(c[ci] < pw.cuts.front() && ci < c.size()) {
295 bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());
297 ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );
301 ret.push_cut(pw.cuts[0]);
302 double prev = pw.cuts[0];
305 while(si < pw.size() && ci <= c.size()) {
306 if(ci == c.size() && prev <= pw.cuts[si]) {
307 ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());
308 ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());
310 }
else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) {
311 if(prev > pw.cuts[si]) {
312 ret.push_seg(portion(pw[si], pw.
segT(prev, si), 1.0));
314 ret.push_seg(pw[si]);
316 ret.push_cut(pw.cuts[si + 1]);
317 prev = pw.cuts[si + 1];
319 }
else if(c[ci] == pw.cuts[si]){
323 ret.
push(elem_portion(pw, si, prev, c[ci]), c[ci]);
330 while(ci < c.size()) {
332 ret.
push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);
350 from = std::min(from, to);
351 to = std::max(temp, to);
353 unsigned i = pw.
segN(from);
355 if(i == pw.size() - 1 || to < pw.cuts[i + 1]) {
356 ret.
push(elem_portion(pw, i, from, to), to);
359 ret.push_seg(portion( pw[i], pw.
segT(from, i), 1.0 ));
361 unsigned fi = pw.
segN(to, i);
363 ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi);
364 ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1);
366 ret.push_seg( portion(pw[fi], 0.0, pw.
segT(to, fi)));
367 if(to != ret.cuts.back()) ret.push_cut(to);
373 Piecewise<T> remove_short_cuts(Piecewise<T>
const &f,
double tol) {
374 if(f.empty())
return f;
376 ret.push_cut(f.cuts[0]);
377 for(
unsigned i=0; i<f.size(); i++){
378 if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {
379 ret.push(f[i], f.cuts[i+1]);
386 Piecewise<T> remove_short_cuts_extending(Piecewise<T>
const &f,
double tol) {
387 if(f.empty())
return f;
389 ret.push_cut(f.cuts[0]);
390 double last = f.cuts[0];
391 for(
unsigned i=0; i<f.size(); i++){
392 if (f.cuts[i+1]-f.cuts[i] >= tol) {
393 ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);
401 std::vector<double> roots(
const Piecewise<T> &pw) {
402 std::vector<double> ret;
403 for(
unsigned i = 0; i < pw.size(); i++) {
404 std::vector<double> sr = roots(pw[i]);
405 for (
unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);
413 Piecewise<T> operator+(Piecewise<T>
const &a,
typename T::output_type b) {
414 boost::function_requires<OffsetableConcept<T> >();
416 Piecewise<T> ret = Piecewise<T>();
418 for(
unsigned i = 0; i < a.size();i++)
419 ret.push_seg(a[i] + b);
423 Piecewise<T> operator-(Piecewise<T>
const &a,
typename T::output_type b) {
424 boost::function_requires<OffsetableConcept<T> >();
426 Piecewise<T> ret = Piecewise<T>();
428 for(
unsigned i = 0; i < a.size();i++)
429 ret.push_seg(a[i] - b);
433 Piecewise<T> operator+=(Piecewise<T>& a,
typename T::output_type b) {
434 boost::function_requires<OffsetableConcept<T> >();
436 if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.);
return a; }
438 for(
unsigned i = 0; i < a.size();i++)
443 Piecewise<T> operator-=(Piecewise<T>& a,
typename T::output_type b) {
444 boost::function_requires<OffsetableConcept<T> >();
446 if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.);
return a; }
448 for(
unsigned i = 0;i < a.size();i++)
455 Piecewise<T> operator-(Piecewise<T>
const &a) {
456 boost::function_requires<ScalableConcept<T> >();
460 for(
unsigned i = 0; i < a.size();i++)
461 ret.push_seg(- a[i]);
465 Piecewise<T> operator*(Piecewise<T>
const &a,
double b) {
466 boost::function_requires<ScalableConcept<T> >();
468 if(a.empty())
return Piecewise<T>();
472 for(
unsigned i = 0; i < a.size();i++)
473 ret.push_seg(a[i] * b);
477 Piecewise<T> operator/(Piecewise<T>
const &a,
double b) {
478 boost::function_requires<ScalableConcept<T> >();
481 if(a.empty())
return Piecewise<T>();
485 for(
unsigned i = 0; i < a.size();i++)
486 ret.push_seg(a[i] / b);
490 Piecewise<T> operator*=(Piecewise<T>& a,
double b) {
491 boost::function_requires<ScalableConcept<T> >();
493 if(a.empty())
return Piecewise<T>();
495 for(
unsigned i = 0; i < a.size();i++)
500 Piecewise<T> operator/=(Piecewise<T>& a,
double b) {
501 boost::function_requires<ScalableConcept<T> >();
504 if(a.empty())
return Piecewise<T>();
506 for(
unsigned i = 0; i < a.size();i++)
513 Piecewise<T> operator+(Piecewise<T>
const &a, Piecewise<T>
const &b) {
514 boost::function_requires<AddableConcept<T> >();
517 Piecewise<T> ret = Piecewise<T>();
518 assert(pa.size() == pb.size());
520 for (
unsigned i = 0; i < pa.size(); i++)
521 ret.push_seg(pa[i] + pb[i]);
525 Piecewise<T> operator-(Piecewise<T>
const &a, Piecewise<T>
const &b) {
526 boost::function_requires<AddableConcept<T> >();
529 Piecewise<T> ret = Piecewise<T>();
530 assert(pa.size() == pb.size());
532 for (
unsigned i = 0; i < pa.size(); i++)
533 ret.push_seg(pa[i] - pb[i]);
537 inline Piecewise<T> operator+=(Piecewise<T> &a, Piecewise<T>
const &b) {
542 inline Piecewise<T> operator-=(Piecewise<T> &a, Piecewise<T>
const &b) {
547 template<
typename T1,
typename T2>
548 Piecewise<T2> operator*(Piecewise<T1>
const &a, Piecewise<T2>
const &b) {
554 Piecewise<T2> ret = Piecewise<T2>();
555 assert(pa.size() == pb.size());
557 for (
unsigned i = 0; i < pa.size(); i++)
558 ret.push_seg(pa[i] * pb[i]);
563 inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T>
const &b) {
568 Piecewise<SBasis> divide(Piecewise<SBasis>
const &a, Piecewise<SBasis>
const &b,
unsigned k);
573 divide(Piecewise<SBasis>
const &a, Piecewise<SBasis>
const &b,
double tol,
unsigned k,
double zero=1.e-3);
575 divide(SBasis
const &a, Piecewise<SBasis>
const &b,
double tol,
unsigned k,
double zero=1.e-3);
577 divide(Piecewise<SBasis>
const &a, SBasis
const &b,
double tol,
unsigned k,
double zero=1.e-3);
579 divide(SBasis
const &a, SBasis
const &b,
double tol,
unsigned k,
double zero=1.e-3);
582 std::map<double,unsigned> compose_pullback(std::vector<double>
const &cuts, SBasis
const &g);
583 int compose_findSegIdx(std::map<double,unsigned>::iterator
const &cut,
584 std::map<double,unsigned>::iterator
const &next,
585 std::vector<double>
const &levels,
590 Piecewise<T> compose(Piecewise<T>
const &f, SBasis
const &g){
592 if (f.empty())
return result;
593 if (g.isZero())
return Piecewise<T>(f(0));
595 double t0 = f.cuts[0], width = f.cuts[1] - t0;
596 return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));
600 Interval bs = bounds_fast(g);
601 if (f.cuts.front() > bs.max() || bs.min() > f.cuts.back()){
602 int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
603 double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
604 return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));
607 std::vector<double> levels;
608 levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);
610 std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);
613 result.cuts.push_back(0.);
614 std::map<double,unsigned>::iterator cut=cuts_pb.begin();
615 std::map<double,unsigned>::iterator next=cut; next++;
616 while(next!=cuts_pb.end()){
620 int idx = compose_findSegIdx(cut,next,levels,g);
621 double t0=(*cut).first;
622 double t1=(*next).first;
624 SBasis sub_g=compose(g, Linear(t0,t1));
625 sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),
626 (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);
627 result.push(compose(f[idx],sub_g),t1);
636 Piecewise<T> compose(Piecewise<T>
const &f, Piecewise<SBasis>
const &g){
638 for(
unsigned i = 0; i < g.segs.size(); i++){
639 Piecewise<T> fgi=compose(f, g.segs[i]);
640 fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));
646 template <
typename T>
647 Piecewise<T> Piecewise<T>::operator()(SBasis f){
return compose((*
this),f);}
648 template <
typename T>
649 Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){
return compose((*
this),f);}
652 Piecewise<T> integral(Piecewise<T>
const &a) {
654 result.segs.resize(a.segs.size());
655 result.cuts = a.cuts;
656 typename T::output_type c = a.segs[0].at0();
657 for(
unsigned i = 0; i < a.segs.size(); i++){
658 result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);
659 result.segs[i]+= c-result.segs[i].at0();
660 c = result.segs[i].at1();
666 Piecewise<T> derivative(Piecewise<T>
const &a) {
668 result.segs.resize(a.segs.size());
669 result.cuts = a.cuts;
670 for(
unsigned i = 0; i < a.segs.size(); i++){
671 result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);
676 std::vector<double> roots(Piecewise<SBasis>
const &f);
680 #endif //SEEN_GEOM_PW_SB_H
Definition: concepts.h:59
void push(const T &s, double to)
Definition: piecewise.h:86
Piecewise< T > partition(const Piecewise< T > &pw, std::vector< double > const &c)
Definition: piecewise.h:276
Definition: piecewise.h:45
unsigned segN(double t, int low=0, int high=-1) const
Definition: piecewise.h:102
double segT(double t, int i=-1) const
Definition: piecewise.h:124