![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: ce1fb12)
raw | patch | inline | side by side (parent: ce1fb12)
| author | joncruz <joncruz@users.sourceforge.net> | |
| Mon, 17 Dec 2007 06:59:40 +0000 (06:59 +0000) | ||
| committer | joncruz <joncruz@users.sourceforge.net> | |
| Mon, 17 Dec 2007 06:59:40 +0000 (06:59 +0000) |
diff --git a/src/2geom/ord.h b/src/2geom/ord.h
index d0e348aece9778f44180d5602494281a1b0bd521..735de6a3ade5cef8d7f4ce4fd263d32e512dba7f 100644 (file)
--- a/src/2geom/ord.h
+++ b/src/2geom/ord.h
namespace {
-enum Cmp {\r
- LESS_THAN=-1,\r
- GREATER_THAN=1,\r
- EQUAL_TO=0\r
-};\r
+enum Cmp {
+ LESS_THAN=-1,
+ GREATER_THAN=1,
+ EQUAL_TO=0
+};
inline Cmp operator-(Cmp x) {
switch(x) {
return EQUAL_TO;
}
}
-\r
-template <typename T1, typename T2>\r
-inline Cmp cmp(T1 const &a, T2 const &b) {\r
- if ( a < b ) {\r
- return LESS_THAN;\r
- } else if ( b < a ) {\r
- return GREATER_THAN;\r
- } else {\r
- return EQUAL_TO;\r
- }\r
+
+template <typename T1, typename T2>
+inline Cmp cmp(T1 const &a, T2 const &b) {
+ if ( a < b ) {
+ return LESS_THAN;
+ } else if ( b < a ) {
+ return GREATER_THAN;
+ } else {
+ return EQUAL_TO;
+ }
}
}
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index c70ecd42cb69be42e8bfd2e4ef846dc190e0ea81..1715997684cec7c7ccff95f8e893c0cd5fe9cdbf 100644 (file)
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
-/*\r
- * piecewise.h - Piecewise function class\r
- *\r
- * Copyright 2007 Michael Sloan <mgsloan@gmail.com>\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, output to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- *\r
- */\r
-\r
-#ifndef SEEN_GEOM_PW_SB_H\r
-#define SEEN_GEOM_PW_SB_H\r
-\r
-#include "sbasis.h"\r
-#include <vector>\r
-#include <map>\r
-\r
-#include "concepts.h"\r
-#include "isnan.h"\r
-#include <boost/concept_check.hpp>\r
-\r
-namespace Geom {\r
-\r
-template <typename T>\r
-class Piecewise {\r
- BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);\r
-\r
- public:\r
- std::vector<double> cuts;\r
- std::vector<T> segs;\r
- //segs[i] stretches from cuts[i] to cuts[i+1].\r
-\r
- Piecewise() {}\r
-\r
- explicit Piecewise(const T &s) {\r
- push_cut(0.);\r
- push_seg(s);\r
- push_cut(1.);\r
- }\r
-\r
- typedef typename T::output_type output_type;\r
-\r
- explicit Piecewise(const output_type & v) {\r
- push_cut(0.);\r
- push_seg(T(v));\r
- push_cut(1.);\r
- }\r
-\r
- inline T operator[](unsigned i) const { return segs[i]; }\r
- inline T &operator[](unsigned i) { return segs[i]; }\r
- inline output_type operator()(double t) const { return valueAt(t); }\r
- inline output_type valueAt(double t) const {\r
- unsigned n = segN(t);\r
- return segs[n](segT(t, n));\r
- }\r
- //TODO: maybe it is not a good idea to have this?\r
- Piecewise<T> operator()(SBasis f);\r
- Piecewise<T> operator()(Piecewise<SBasis>f);\r
-\r
- inline unsigned size() const { return segs.size(); }\r
- inline bool empty() const { return segs.empty(); }\r
-\r
- /**Convenience/implementation hiding function to add segment/cut pairs.\r
- * Asserts that basic size and order invariants are correct\r
- */\r
- inline void push(const T &s, double to) {\r
- assert(cuts.size() - segs.size() == 1);\r
- push_seg(s);\r
- push_cut(to);\r
- }\r
- //Convenience/implementation hiding function to add cuts.\r
- inline void push_cut(double c) {\r
- assert_invariants(cuts.empty() || c > cuts.back()); \r
- cuts.push_back(c);\r
- }\r
- //Convenience/implementation hiding function to add segments.\r
- inline void push_seg(const T &s) { segs.push_back(s); }\r
-\r
- /**Returns the segment index which corresponds to a 'global' piecewise time.\r
- * Also takes optional low/high parameters to expedite the search for the segment.\r
- */\r
- inline unsigned segN(double t, int low = 0, int high = -1) const {\r
- high = (high == -1) ? size() : high;\r
- if(t < cuts[0]) return 0;\r
- if(t >= cuts[size()]) return size() - 1;\r
- while(low < high) {\r
- int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences\r
- double mv = cuts[mid];\r
- if(mv < t) {\r
- if(t < cuts[mid + 1]) return mid; else low = mid + 1;\r
- } else if(t < mv) {\r
- if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;\r
- } else {\r
- return mid;\r
- }\r
- }\r
- return low;\r
- }\r
-\r
- /**Returns the time within a segment, given the 'global' piecewise time.\r
- * Also takes an optional index parameter which may be used for efficiency or to find the time on a\r
- * segment outside its range. If it is left to its default, -1, it will call segN to find the index.\r
- */\r
- inline double segT(double t, int i = -1) const {\r
- if(i == -1) i = segN(t);\r
- assert(i >= 0);\r
- return (t - cuts[i]) / (cuts[i+1] - cuts[i]);\r
- }\r
-\r
- inline double mapToDomain(double t, unsigned i) const {\r
- return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]\r
- } \r
-\r
- //Offsets the piecewise domain\r
- inline void offsetDomain(double o) {\r
- assert(is_finite(o));\r
- if(o != 0)\r
- for(unsigned i = 0; i <= size(); i++)\r
- cuts[i] += o;\r
- }\r
-\r
- //Scales the domain of the function by a value. 0 will result in an empty Piecewise.\r
- inline void scaleDomain(double s) {\r
- assert(s > 0);\r
- if(s == 0) {\r
- cuts.clear(); segs.clear();\r
- return;\r
- }\r
- for(unsigned i = 0; i <= size(); i++)\r
- cuts[i] *= s;\r
- }\r
-\r
- //Retrieves the domain in interval form\r
- inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }\r
-\r
- //Transforms the domain into another interval\r
- inline void setDomain(Interval dom) {\r
- if(empty()) return;\r
- if(dom.isEmpty()) {\r
- cuts.clear(); segs.clear();\r
- return;\r
- }\r
- double cf = cuts.front();\r
- double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);\r
- for(unsigned i = 0; i <= size(); i++)\r
- cuts[i] = (cuts[i] - cf) * s + o;\r
- }\r
-\r
- //Concatenates this Piecewise function with another, offseting time of the other to match the end.\r
- inline void concat(const Piecewise<T> &other) {\r
- if(other.empty()) return;\r
-\r
- if(empty()) { \r
- cuts = other.cuts; segs = other.segs;\r
- return;\r
- }\r
-\r
- segs.insert(segs.end(), other.segs.begin(), other.segs.end());\r
- double t = cuts.back() - other.cuts.front(); \r
- for(unsigned i = 0; i < other.size(); i++)\r
- push_cut(other.cuts[i + 1] + t);\r
- }\r
-\r
- //Like concat, but ensures continuity.\r
- inline void continuousConcat(const Piecewise<T> &other) {\r
- boost::function_requires<AddableConcept<typename T::output_type> >();\r
- if(other.empty()) return;\r
- typename T::output_type y = segs.back().at1() - other.segs.front().at0();\r
-\r
- if(empty()) {\r
- for(unsigned i = 0; i < other.size(); i++)\r
- push_seg(other[i] + y);\r
- cuts = other.cuts;\r
- return;\r
- }\r
-\r
- double t = cuts.back() - other.cuts.front();\r
- for(unsigned i = 0; i < other.size(); i++)\r
- push(other[i] + y, other.cuts[i + 1] + t);\r
- }\r
-\r
- //returns true if the Piecewise<T> meets some basic invariants.\r
- inline bool invariants() const {\r
- // segs between cuts\r
- if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))\r
- return false;\r
- // cuts in order\r
- for(unsigned i = 0; i < segs.size(); i++)\r
- if(cuts[i] >= cuts[i+1])\r
- return false;\r
- return true;\r
- }\r
-\r
-};\r
-\r
-template<typename T>\r
-inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {\r
- boost::function_requires<FragmentConcept<T> >();\r
- \r
- if(f.empty()) return typename FragmentConcept<T>::BoundsType();\r
- typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));\r
- for(unsigned i = 1; i < f.size(); i++)\r
- ret.unionWith(bounds_fast(f[i])); \r
- return ret;\r
-}\r
-\r
-template<typename T>\r
-inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {\r
- boost::function_requires<FragmentConcept<T> >();\r
- \r
- if(f.empty()) return typename FragmentConcept<T>::BoundsType();\r
- typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));\r
- for(unsigned i = 1; i < f.size(); i++)\r
- ret.unionWith(bounds_exact(f[i])); \r
- return ret;\r
-}\r
-\r
-template<typename T>\r
-inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const Interval &m) {\r
- boost::function_requires<FragmentConcept<T> >();\r
- \r
- if(f.empty()) return typename FragmentConcept<T>::BoundsType();\r
- if(m.isEmpty()) return typename FragmentConcept<T>::BoundsType(f(m.min()));\r
- \r
- unsigned fi = f.segN(m.min()), ti = f.segN(m.max());\r
- double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);\r
- \r
- if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));\r
- \r
- typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));\r
- for(unsigned i = fi + 1; i < ti; i++)\r
- ret.unionWith(bounds_exact(f[i]));\r
- if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));\r
- \r
- return ret;\r
-}\r
-\r
-//returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.\r
-template<typename T>\r
-T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {\r
- assert(i < a.size());\r
- double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);\r
- return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );\r
-}\r
-\r
-/**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);\r
- * Further subdivides the Piecewise<T> such that there is a cut at every value in c.\r
- * Precondition: c sorted lower to higher.\r
- * \r
- * //Given Piecewise<T> a and b:\r
- * Piecewise<T> ac = a.partition(b.cuts);\r
- * Piecewise<T> bc = b.partition(a.cuts);\r
- * //ac.cuts should be equivalent to bc.cuts\r
- */\r
-template<typename T>\r
-Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {\r
- assert(pw.invariants());\r
- if(c.empty()) return Piecewise<T>(pw);\r
-\r
- Piecewise<T> ret = Piecewise<T>();\r
- ret.cuts.reserve(c.size() + pw.cuts.size());\r
- ret.segs.reserve(c.size() + pw.cuts.size() - 1);\r
-\r
- if(pw.empty()) {\r
- ret.cuts = c;\r
- for(unsigned i = 0; i < c.size() - 1; i++)\r
- ret.push_seg(T());\r
- return ret;\r
- }\r
-\r
- unsigned si = 0, ci = 0; //Segment index, Cut index\r
-\r
- //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment\r
- while(c[ci] < pw.cuts.front() && ci < c.size()) {\r
- bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());\r
- ret.push_cut(c[ci]);\r
- ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );\r
- ci++;\r
- }\r
-\r
- ret.push_cut(pw.cuts[0]);\r
- double prev = pw.cuts[0]; //previous cut\r
- //Loop which handles cuts within the Piecewise<T> domain\r
- //Should have the cuts = segs + 1 invariant\r
- while(si < pw.size() && ci <= c.size()) {\r
- if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest\r
- ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());\r
- ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());\r
- return ret;\r
- }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) { //no more cuts within this segment, finalize\r
- if(prev > pw.cuts[si]) { //segment already has cuts, so portion is required\r
- ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));\r
- } else { //plain copy is fine\r
- ret.push_seg(pw[si]);\r
- }\r
- ret.push_cut(pw.cuts[si + 1]);\r
- prev = pw.cuts[si + 1];\r
- si++;\r
- } else if(c[ci] == pw.cuts[si]){ //coincident\r
- //Already finalized the seg with the code immediately above\r
- ci++;\r
- } else { //plain old subdivision\r
- ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);\r
- prev = c[ci];\r
- ci++;\r
- }\r
- }\r
- \r
- //input cuts extend further than this Piecewise<T>, extend the last segment.\r
- while(ci < c.size()) {\r
- if(c[ci] > prev) {\r
- ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);\r
- prev = c[ci];\r
- }\r
- ci++;\r
- }\r
- return ret;\r
-}\r
-\r
-/**Piecewise<T> portion(const Piecewise<T> &pw, double from, double to);\r
- * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].\r
- */\r
-template<typename T>\r
-Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {\r
- if(pw.empty() || from == to) return Piecewise<T>();\r
-\r
- Piecewise<T> ret;\r
-\r
- double temp = from;\r
- from = std::min(from, to);\r
- to = std::max(temp, to);\r
- \r
- unsigned i = pw.segN(from);\r
- ret.push_cut(from);\r
- if(i == pw.size() - 1 || to < pw.cuts[i + 1]) { //to/from inhabit the same segment\r
- ret.push(elem_portion(pw, i, from, to), to);\r
- return ret;\r
- }\r
- ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));\r
- i++;\r
- unsigned fi = pw.segN(to, i);\r
-\r
- ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi); //copy segs\r
- ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1); //and their cuts\r
-\r
- ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));\r
- if(to != ret.cuts.back()) ret.push_cut(to);\r
- ret.invariants();\r
- return ret;\r
-}\r
-\r
-template<typename T>\r
-Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {\r
- if(f.empty()) return f;\r
- Piecewise<T> ret;\r
- ret.push_cut(f.cuts[0]);\r
- for(unsigned i=0; i<f.size(); i++){\r
- if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {\r
- ret.push(f[i], f.cuts[i+1]);\r
- }\r
- }\r
- return ret;\r
-}\r
-\r
-template<typename T>\r
-Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {\r
- if(f.empty()) return f;\r
- Piecewise<T> ret;\r
- ret.push_cut(f.cuts[0]);\r
- double last = f.cuts[0]; // last cut included\r
- for(unsigned i=0; i<f.size(); i++){\r
- if (f.cuts[i+1]-f.cuts[i] >= tol) {\r
- ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);\r
- last = f.cuts[i+1];\r
- }\r
- }\r
- return ret;\r
-}\r
-\r
-template<typename T>\r
-std::vector<double> roots(const Piecewise<T> &pw) {\r
- std::vector<double> ret;\r
- for(unsigned i = 0; i < pw.size(); i++) {\r
- std::vector<double> sr = roots(pw[i]);\r
- for (unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);\r
-\r
- }\r
- return ret;\r
-}\r
-\r
-//IMPL: OffsetableConcept\r
-template<typename T>\r
-Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {\r
- boost::function_requires<OffsetableConcept<T> >();\r
-//TODO:empty\r
- Piecewise<T> ret = Piecewise<T>();\r
- ret.cuts = a.cuts;\r
- for(unsigned i = 0; i < a.size();i++)\r
- ret.push_seg(a[i] + b);\r
- return ret;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {\r
- boost::function_requires<OffsetableConcept<T> >();\r
-//TODO: empty\r
- Piecewise<T> ret = Piecewise<T>();\r
- ret.cuts = a.cuts;\r
- for(unsigned i = 0; i < a.size();i++)\r
- ret.push_seg(a[i] - b);\r
- return ret;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator+=(Piecewise<T>& a, typename T::output_type b) {\r
- boost::function_requires<OffsetableConcept<T> >();\r
-\r
- if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }\r
-\r
- for(unsigned i = 0; i < a.size();i++)\r
- a[i] += b;\r
- return a;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator-=(Piecewise<T>& a, typename T::output_type b) {\r
- boost::function_requires<OffsetableConcept<T> >();\r
-\r
- if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }\r
-\r
- for(unsigned i = 0;i < a.size();i++)\r
- a[i] -= b;\r
- return a;\r
-}\r
-\r
-//IMPL: ScalableConcept\r
-template<typename T>\r
-Piecewise<T> operator-(Piecewise<T> const &a) {\r
- boost::function_requires<ScalableConcept<T> >();\r
-\r
- Piecewise<T> ret;\r
- ret.cuts = a.cuts;\r
- for(unsigned i = 0; i < a.size();i++)\r
- ret.push_seg(- a[i]);\r
- return ret;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator*(Piecewise<T> const &a, double b) {\r
- boost::function_requires<ScalableConcept<T> >();\r
-\r
- if(a.empty()) return Piecewise<T>();\r
-\r
- Piecewise<T> ret;\r
- ret.cuts = a.cuts;\r
- for(unsigned i = 0; i < a.size();i++)\r
- ret.push_seg(a[i] * b);\r
- return ret;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator/(Piecewise<T> const &a, double b) {\r
- boost::function_requires<ScalableConcept<T> >();\r
-\r
- //FIXME: b == 0?\r
- if(a.empty()) return Piecewise<T>();\r
-\r
- Piecewise<T> ret;\r
- ret.cuts = a.cuts;\r
- for(unsigned i = 0; i < a.size();i++)\r
- ret.push_seg(a[i] / b);\r
- return ret;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator*=(Piecewise<T>& a, double b) {\r
- boost::function_requires<ScalableConcept<T> >();\r
-\r
- if(a.empty()) return Piecewise<T>();\r
-\r
- for(unsigned i = 0; i < a.size();i++)\r
- a[i] *= b;\r
- return a;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator/=(Piecewise<T>& a, double b) {\r
- boost::function_requires<ScalableConcept<T> >();\r
-\r
- //FIXME: b == 0?\r
- if(a.empty()) return Piecewise<T>();\r
-\r
- for(unsigned i = 0; i < a.size();i++)\r
- a[i] /= b;\r
- return a;\r
-}\r
-\r
-//IMPL: AddableConcept\r
-template<typename T>\r
-Piecewise<T> operator+(Piecewise<T> const &a, Piecewise<T> const &b) {\r
- boost::function_requires<AddableConcept<T> >();\r
-\r
- Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);\r
- Piecewise<T> ret = Piecewise<T>();\r
- assert(pa.size() == pb.size());\r
- ret.cuts = pa.cuts;\r
- for (unsigned i = 0; i < pa.size(); i++)\r
- ret.push_seg(pa[i] + pb[i]);\r
- return ret;\r
-}\r
-template<typename T>\r
-Piecewise<T> operator-(Piecewise<T> const &a, Piecewise<T> const &b) {\r
- boost::function_requires<AddableConcept<T> >();\r
-\r
- Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);\r
- Piecewise<T> ret = Piecewise<T>();\r
- assert(pa.size() == pb.size());\r
- ret.cuts = pa.cuts;\r
- for (unsigned i = 0; i < pa.size(); i++)\r
- ret.push_seg(pa[i] - pb[i]);\r
- return ret;\r
-}\r
-template<typename T>\r
-inline Piecewise<T> operator+=(Piecewise<T> &a, Piecewise<T> const &b) {\r
- a = a+b;\r
- return a;\r
-}\r
-template<typename T>\r
-inline Piecewise<T> operator-=(Piecewise<T> &a, Piecewise<T> const &b) {\r
- a = a-b;\r
- return a;\r
-}\r
-\r
-template<typename T1,typename T2>\r
-Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {\r
- //function_requires<MultiplicableConcept<T1> >();\r
- //function_requires<MultiplicableConcept<T2> >();\r
-\r
- Piecewise<T1> pa = partition(a, b.cuts);\r
- Piecewise<T2> pb = partition(b, a.cuts);\r
- Piecewise<T2> ret = Piecewise<T2>();\r
- assert(pa.size() == pb.size());\r
- ret.cuts = pa.cuts;\r
- for (unsigned i = 0; i < pa.size(); i++)\r
- ret.push_seg(pa[i] * pb[i]);\r
- return ret;\r
-}\r
-\r
-template<typename T>\r
-inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T> const &b) {\r
- a = a * b;\r
- return a;\r
-}\r
-\r
-Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);\r
-//TODO: replace divide(a,b,k) by divide(a,b,tol,k)?\r
-//TODO: atm, relative error is ~(tol/a)%. Find a way to make it independant of a.\r
-//Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero). \r
-Piecewise<SBasis> \r
-divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);\r
-Piecewise<SBasis> \r
-divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);\r
-Piecewise<SBasis> \r
-divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);\r
-Piecewise<SBasis> \r
-divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);\r
-\r
-//Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.\r
-std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);\r
-int compose_findSegIdx(std::map<double,unsigned>::iterator const &cut,\r
- std::map<double,unsigned>::iterator const &next,\r
- std::vector<double> const &levels,\r
- SBasis const &g);\r
-\r
-//TODO: add concept check\r
-template<typename T>\r
-Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){\r
- Piecewise<T> result;\r
- if (f.empty()) return result;\r
- if (g.isZero()) return Piecewise<T>(f(0));\r
- if (f.size()==1){\r
- double t0 = f.cuts[0], width = f.cuts[1] - t0;\r
- return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));\r
- }\r
- \r
- //first check bounds...\r
- Interval bs = bounds_fast(g);\r
- if (f.cuts.front() > bs.max() || bs.min() > f.cuts.back()){\r
- int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;\r
- double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;\r
- return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g)); \r
- }\r
- \r
- std::vector<double> levels;//we can forget first and last cuts...\r
- levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);\r
- //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.\r
- std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);\r
- \r
- //-- Compose each piece of g with the relevant seg of f.\r
- result.cuts.push_back(0.);\r
- std::map<double,unsigned>::iterator cut=cuts_pb.begin();\r
- std::map<double,unsigned>::iterator next=cut; next++;\r
- while(next!=cuts_pb.end()){\r
- //assert(std::abs(int((*cut).second-(*next).second))<1);\r
- //TODO: find a way to recover from this error? the root finder missed some root;\r
- // the levels/variations of f might be too close/fast...\r
- int idx = compose_findSegIdx(cut,next,levels,g);\r
- double t0=(*cut).first;\r
- double t1=(*next).first;\r
- \r
- SBasis sub_g=compose(g, Linear(t0,t1));\r
- sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),\r
- (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);\r
- result.push(compose(f[idx],sub_g),t1);\r
- cut++;\r
- next++;\r
- }\r
- return(result);\r
-} \r
-\r
-//TODO: add concept check for following composition functions\r
-template<typename T>\r
-Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){\r
- Piecewise<T> result;\r
- for(unsigned i = 0; i < g.segs.size(); i++){\r
- Piecewise<T> fgi=compose(f, g.segs[i]);\r
- fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));\r
- result.concat(fgi);\r
- }\r
- return result;\r
-}\r
-\r
-template <typename T>\r
-Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}\r
-template <typename T>\r
-Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}\r
-\r
-template<typename T>\r
-Piecewise<T> integral(Piecewise<T> const &a) {\r
- Piecewise<T> result;\r
- result.segs.resize(a.segs.size());\r
- result.cuts = a.cuts;\r
- typename T::output_type c = a.segs[0].at0();\r
- for(unsigned i = 0; i < a.segs.size(); i++){\r
- result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);\r
- result.segs[i]+= c-result.segs[i].at0();\r
- c = result.segs[i].at1();\r
- }\r
- return result;\r
-}\r
-\r
-template<typename T>\r
-Piecewise<T> derivative(Piecewise<T> const &a) {\r
- Piecewise<T> result;\r
- result.segs.resize(a.segs.size());\r
- result.cuts = a.cuts;\r
- for(unsigned i = 0; i < a.segs.size(); i++){\r
- result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);\r
- }\r
- return result;\r
-}\r
-\r
-std::vector<double> roots(Piecewise<SBasis> const &f);\r
-\r
-}\r
-\r
-#endif //SEEN_GEOM_PW_SB_H\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * piecewise.h - Piecewise function class
+ *
+ * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, output to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ */
+
+#ifndef SEEN_GEOM_PW_SB_H
+#define SEEN_GEOM_PW_SB_H
+
+#include "sbasis.h"
+#include <vector>
+#include <map>
+
+#include "concepts.h"
+#include "isnan.h"
+#include <boost/concept_check.hpp>
+
+namespace Geom {
+
+template <typename T>
+class Piecewise {
+ BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);
+
+ public:
+ std::vector<double> cuts;
+ std::vector<T> segs;
+ //segs[i] stretches from cuts[i] to cuts[i+1].
+
+ Piecewise() {}
+
+ explicit Piecewise(const T &s) {
+ push_cut(0.);
+ push_seg(s);
+ push_cut(1.);
+ }
+
+ typedef typename T::output_type output_type;
+
+ explicit Piecewise(const output_type & v) {
+ push_cut(0.);
+ push_seg(T(v));
+ push_cut(1.);
+ }
+
+ inline T operator[](unsigned i) const { return segs[i]; }
+ inline T &operator[](unsigned i) { return segs[i]; }
+ inline output_type operator()(double t) const { return valueAt(t); }
+ inline output_type valueAt(double t) const {
+ unsigned n = segN(t);
+ return segs[n](segT(t, n));
+ }
+ //TODO: maybe it is not a good idea to have this?
+ Piecewise<T> operator()(SBasis f);
+ Piecewise<T> operator()(Piecewise<SBasis>f);
+
+ inline unsigned size() const { return segs.size(); }
+ inline bool empty() const { return segs.empty(); }
+
+ /**Convenience/implementation hiding function to add segment/cut pairs.
+ * Asserts that basic size and order invariants are correct
+ */
+ inline void push(const T &s, double to) {
+ assert(cuts.size() - segs.size() == 1);
+ push_seg(s);
+ push_cut(to);
+ }
+ //Convenience/implementation hiding function to add cuts.
+ inline void push_cut(double c) {
+ assert_invariants(cuts.empty() || c > cuts.back());
+ cuts.push_back(c);
+ }
+ //Convenience/implementation hiding function to add segments.
+ inline void push_seg(const T &s) { segs.push_back(s); }
+
+ /**Returns the segment index which corresponds to a 'global' piecewise time.
+ * Also takes optional low/high parameters to expedite the search for the segment.
+ */
+ inline unsigned segN(double t, int low = 0, int high = -1) const {
+ high = (high == -1) ? size() : high;
+ if(t < cuts[0]) return 0;
+ if(t >= cuts[size()]) return size() - 1;
+ while(low < high) {
+ int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences
+ double mv = cuts[mid];
+ if(mv < t) {
+ if(t < cuts[mid + 1]) return mid; else low = mid + 1;
+ } else if(t < mv) {
+ if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;
+ } else {
+ return mid;
+ }
+ }
+ return low;
+ }
+
+ /**Returns the time within a segment, given the 'global' piecewise time.
+ * Also takes an optional index parameter which may be used for efficiency or to find the time on a
+ * segment outside its range. If it is left to its default, -1, it will call segN to find the index.
+ */
+ inline double segT(double t, int i = -1) const {
+ if(i == -1) i = segN(t);
+ assert(i >= 0);
+ return (t - cuts[i]) / (cuts[i+1] - cuts[i]);
+ }
+
+ inline double mapToDomain(double t, unsigned i) const {
+ return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]
+ }
+
+ //Offsets the piecewise domain
+ inline void offsetDomain(double o) {
+ assert(is_finite(o));
+ if(o != 0)
+ for(unsigned i = 0; i <= size(); i++)
+ cuts[i] += o;
+ }
+
+ //Scales the domain of the function by a value. 0 will result in an empty Piecewise.
+ inline void scaleDomain(double s) {
+ assert(s > 0);
+ if(s == 0) {
+ cuts.clear(); segs.clear();
+ return;
+ }
+ for(unsigned i = 0; i <= size(); i++)
+ cuts[i] *= s;
+ }
+
+ //Retrieves the domain in interval form
+ inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }
+
+ //Transforms the domain into another interval
+ inline void setDomain(Interval dom) {
+ if(empty()) return;
+ if(dom.isEmpty()) {
+ cuts.clear(); segs.clear();
+ return;
+ }
+ double cf = cuts.front();
+ double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
+ for(unsigned i = 0; i <= size(); i++)
+ cuts[i] = (cuts[i] - cf) * s + o;
+ }
+
+ //Concatenates this Piecewise function with another, offseting time of the other to match the end.
+ inline void concat(const Piecewise<T> &other) {
+ if(other.empty()) return;
+
+ if(empty()) {
+ cuts = other.cuts; segs = other.segs;
+ return;
+ }
+
+ segs.insert(segs.end(), other.segs.begin(), other.segs.end());
+ double t = cuts.back() - other.cuts.front();
+ for(unsigned i = 0; i < other.size(); i++)
+ push_cut(other.cuts[i + 1] + t);
+ }
+
+ //Like concat, but ensures continuity.
+ inline void continuousConcat(const Piecewise<T> &other) {
+ boost::function_requires<AddableConcept<typename T::output_type> >();
+ if(other.empty()) return;
+ typename T::output_type y = segs.back().at1() - other.segs.front().at0();
+
+ if(empty()) {
+ for(unsigned i = 0; i < other.size(); i++)
+ push_seg(other[i] + y);
+ cuts = other.cuts;
+ return;
+ }
+
+ double t = cuts.back() - other.cuts.front();
+ for(unsigned i = 0; i < other.size(); i++)
+ push(other[i] + y, other.cuts[i + 1] + t);
+ }
+
+ //returns true if the Piecewise<T> meets some basic invariants.
+ inline bool invariants() const {
+ // segs between cuts
+ if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))
+ return false;
+ // cuts in order
+ for(unsigned i = 0; i < segs.size(); i++)
+ if(cuts[i] >= cuts[i+1])
+ return false;
+ return true;
+ }
+
+};
+
+template<typename T>
+inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {
+ boost::function_requires<FragmentConcept<T> >();
+
+ if(f.empty()) return typename FragmentConcept<T>::BoundsType();
+ typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));
+ for(unsigned i = 1; i < f.size(); i++)
+ ret.unionWith(bounds_fast(f[i]));
+ return ret;
+}
+
+template<typename T>
+inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {
+ boost::function_requires<FragmentConcept<T> >();
+
+ if(f.empty()) return typename FragmentConcept<T>::BoundsType();
+ typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));
+ for(unsigned i = 1; i < f.size(); i++)
+ ret.unionWith(bounds_exact(f[i]));
+ return ret;
+}
+
+template<typename T>
+inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const Interval &m) {
+ boost::function_requires<FragmentConcept<T> >();
+
+ if(f.empty()) return typename FragmentConcept<T>::BoundsType();
+ if(m.isEmpty()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
+
+ unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
+ double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
+
+ if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));
+
+ typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));
+ for(unsigned i = fi + 1; i < ti; i++)
+ ret.unionWith(bounds_exact(f[i]));
+ if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));
+
+ return ret;
+}
+
+//returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.
+template<typename T>
+T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {
+ assert(i < a.size());
+ double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);
+ return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );
+}
+
+/**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);
+ * Further subdivides the Piecewise<T> such that there is a cut at every value in c.
+ * Precondition: c sorted lower to higher.
+ *
+ * //Given Piecewise<T> a and b:
+ * Piecewise<T> ac = a.partition(b.cuts);
+ * Piecewise<T> bc = b.partition(a.cuts);
+ * //ac.cuts should be equivalent to bc.cuts
+ */
+template<typename T>
+Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {
+ assert(pw.invariants());
+ if(c.empty()) return Piecewise<T>(pw);
+
+ Piecewise<T> ret = Piecewise<T>();
+ ret.cuts.reserve(c.size() + pw.cuts.size());
+ ret.segs.reserve(c.size() + pw.cuts.size() - 1);
+
+ if(pw.empty()) {
+ ret.cuts = c;
+ for(unsigned i = 0; i < c.size() - 1; i++)
+ ret.push_seg(T());
+ return ret;
+ }
+
+ unsigned si = 0, ci = 0; //Segment index, Cut index
+
+ //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment
+ while(c[ci] < pw.cuts.front() && ci < c.size()) {
+ bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());
+ ret.push_cut(c[ci]);
+ ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );
+ ci++;
+ }
+
+ ret.push_cut(pw.cuts[0]);
+ double prev = pw.cuts[0]; //previous cut
+ //Loop which handles cuts within the Piecewise<T> domain
+ //Should have the cuts = segs + 1 invariant
+ while(si < pw.size() && ci <= c.size()) {
+ if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest
+ ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());
+ ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());
+ return ret;
+ }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) { //no more cuts within this segment, finalize
+ if(prev > pw.cuts[si]) { //segment already has cuts, so portion is required
+ ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));
+ } else { //plain copy is fine
+ ret.push_seg(pw[si]);
+ }
+ ret.push_cut(pw.cuts[si + 1]);
+ prev = pw.cuts[si + 1];
+ si++;
+ } else if(c[ci] == pw.cuts[si]){ //coincident
+ //Already finalized the seg with the code immediately above
+ ci++;
+ } else { //plain old subdivision
+ ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);
+ prev = c[ci];
+ ci++;
+ }
+ }
+
+ //input cuts extend further than this Piecewise<T>, extend the last segment.
+ while(ci < c.size()) {
+ if(c[ci] > prev) {
+ ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);
+ prev = c[ci];
+ }
+ ci++;
+ }
+ return ret;
+}
+
+/**Piecewise<T> portion(const Piecewise<T> &pw, double from, double to);
+ * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].
+ */
+template<typename T>
+Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {
+ if(pw.empty() || from == to) return Piecewise<T>();
+
+ Piecewise<T> ret;
+
+ double temp = from;
+ from = std::min(from, to);
+ to = std::max(temp, to);
+
+ unsigned i = pw.segN(from);
+ ret.push_cut(from);
+ if(i == pw.size() - 1 || to < pw.cuts[i + 1]) { //to/from inhabit the same segment
+ ret.push(elem_portion(pw, i, from, to), to);
+ return ret;
+ }
+ ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));
+ i++;
+ unsigned fi = pw.segN(to, i);
+
+ ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi); //copy segs
+ ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1); //and their cuts
+
+ ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));
+ if(to != ret.cuts.back()) ret.push_cut(to);
+ ret.invariants();
+ return ret;
+}
+
+template<typename T>
+Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {
+ if(f.empty()) return f;
+ Piecewise<T> ret;
+ ret.push_cut(f.cuts[0]);
+ for(unsigned i=0; i<f.size(); i++){
+ if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {
+ ret.push(f[i], f.cuts[i+1]);
+ }
+ }
+ return ret;
+}
+
+template<typename T>
+Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {
+ if(f.empty()) return f;
+ Piecewise<T> ret;
+ ret.push_cut(f.cuts[0]);
+ double last = f.cuts[0]; // last cut included
+ for(unsigned i=0; i<f.size(); i++){
+ if (f.cuts[i+1]-f.cuts[i] >= tol) {
+ ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);
+ last = f.cuts[i+1];
+ }
+ }
+ return ret;
+}
+
+template<typename T>
+std::vector<double> roots(const Piecewise<T> &pw) {
+ std::vector<double> ret;
+ for(unsigned i = 0; i < pw.size(); i++) {
+ std::vector<double> sr = roots(pw[i]);
+ for (unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);
+
+ }
+ return ret;
+}
+
+//IMPL: OffsetableConcept
+template<typename T>
+Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {
+ boost::function_requires<OffsetableConcept<T> >();
+//TODO:empty
+ Piecewise<T> ret = Piecewise<T>();
+ ret.cuts = a.cuts;
+ for(unsigned i = 0; i < a.size();i++)
+ ret.push_seg(a[i] + b);
+ return ret;
+}
+template<typename T>
+Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {
+ boost::function_requires<OffsetableConcept<T> >();
+//TODO: empty
+ Piecewise<T> ret = Piecewise<T>();
+ ret.cuts = a.cuts;
+ for(unsigned i = 0; i < a.size();i++)
+ ret.push_seg(a[i] - b);
+ return ret;
+}
+template<typename T>
+Piecewise<T> operator+=(Piecewise<T>& a, typename T::output_type b) {
+ boost::function_requires<OffsetableConcept<T> >();
+
+ if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
+
+ for(unsigned i = 0; i < a.size();i++)
+ a[i] += b;
+ return a;
+}
+template<typename T>
+Piecewise<T> operator-=(Piecewise<T>& a, typename T::output_type b) {
+ boost::function_requires<OffsetableConcept<T> >();
+
+ if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
+
+ for(unsigned i = 0;i < a.size();i++)
+ a[i] -= b;
+ return a;
+}
+
+//IMPL: ScalableConcept
+template<typename T>
+Piecewise<T> operator-(Piecewise<T> const &a) {
+ boost::function_requires<ScalableConcept<T> >();
+
+ Piecewise<T> ret;
+ ret.cuts = a.cuts;
+ for(unsigned i = 0; i < a.size();i++)
+ ret.push_seg(- a[i]);
+ return ret;
+}
+template<typename T>
+Piecewise<T> operator*(Piecewise<T> const &a, double b) {
+ boost::function_requires<ScalableConcept<T> >();
+
+ if(a.empty()) return Piecewise<T>();
+
+ Piecewise<T> ret;
+ ret.cuts = a.cuts;
+ for(unsigned i = 0; i < a.size();i++)
+ ret.push_seg(a[i] * b);
+ return ret;
+}
+template<typename T>
+Piecewise<T> operator/(Piecewise<T> const &a, double b) {
+ boost::function_requires<ScalableConcept<T> >();
+
+ //FIXME: b == 0?
+ if(a.empty()) return Piecewise<T>();
+
+ Piecewise<T> ret;
+ ret.cuts = a.cuts;
+ for(unsigned i = 0; i < a.size();i++)
+ ret.push_seg(a[i] / b);
+ return ret;
+}
+template<typename T>
+Piecewise<T> operator*=(Piecewise<T>& a, double b) {
+ boost::function_requires<ScalableConcept<T> >();
+
+ if(a.empty()) return Piecewise<T>();
+
+ for(unsigned i = 0; i < a.size();i++)
+ a[i] *= b;
+ return a;
+}
+template<typename T>
+Piecewise<T> operator/=(Piecewise<T>& a, double b) {
+ boost::function_requires<ScalableConcept<T> >();
+
+ //FIXME: b == 0?
+ if(a.empty()) return Piecewise<T>();
+
+ for(unsigned i = 0; i < a.size();i++)
+ a[i] /= b;
+ return a;
+}
+
+//IMPL: AddableConcept
+template<typename T>
+Piecewise<T> operator+(Piecewise<T> const &a, Piecewise<T> const &b) {
+ boost::function_requires<AddableConcept<T> >();
+
+ Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
+ Piecewise<T> ret = Piecewise<T>();
+ assert(pa.size() == pb.size());
+ ret.cuts = pa.cuts;
+ for (unsigned i = 0; i < pa.size(); i++)
+ ret.push_seg(pa[i] + pb[i]);
+ return ret;
+}
+template<typename T>
+Piecewise<T> operator-(Piecewise<T> const &a, Piecewise<T> const &b) {
+ boost::function_requires<AddableConcept<T> >();
+
+ Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
+ Piecewise<T> ret = Piecewise<T>();
+ assert(pa.size() == pb.size());
+ ret.cuts = pa.cuts;
+ for (unsigned i = 0; i < pa.size(); i++)
+ ret.push_seg(pa[i] - pb[i]);
+ return ret;
+}
+template<typename T>
+inline Piecewise<T> operator+=(Piecewise<T> &a, Piecewise<T> const &b) {
+ a = a+b;
+ return a;
+}
+template<typename T>
+inline Piecewise<T> operator-=(Piecewise<T> &a, Piecewise<T> const &b) {
+ a = a-b;
+ return a;
+}
+
+template<typename T1,typename T2>
+Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {
+ //function_requires<MultiplicableConcept<T1> >();
+ //function_requires<MultiplicableConcept<T2> >();
+
+ Piecewise<T1> pa = partition(a, b.cuts);
+ Piecewise<T2> pb = partition(b, a.cuts);
+ Piecewise<T2> ret = Piecewise<T2>();
+ assert(pa.size() == pb.size());
+ ret.cuts = pa.cuts;
+ for (unsigned i = 0; i < pa.size(); i++)
+ ret.push_seg(pa[i] * pb[i]);
+ return ret;
+}
+
+template<typename T>
+inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T> const &b) {
+ a = a * b;
+ return a;
+}
+
+Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);
+//TODO: replace divide(a,b,k) by divide(a,b,tol,k)?
+//TODO: atm, relative error is ~(tol/a)%. Find a way to make it independant of a.
+//Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero).
+Piecewise<SBasis>
+divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
+Piecewise<SBasis>
+divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
+Piecewise<SBasis>
+divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
+Piecewise<SBasis>
+divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
+
+//Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.
+std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);
+int compose_findSegIdx(std::map<double,unsigned>::iterator const &cut,
+ std::map<double,unsigned>::iterator const &next,
+ std::vector<double> const &levels,
+ SBasis const &g);
+
+//TODO: add concept check
+template<typename T>
+Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){
+ Piecewise<T> result;
+ if (f.empty()) return result;
+ if (g.isZero()) return Piecewise<T>(f(0));
+ if (f.size()==1){
+ double t0 = f.cuts[0], width = f.cuts[1] - t0;
+ return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));
+ }
+
+ //first check bounds...
+ Interval bs = bounds_fast(g);
+ if (f.cuts.front() > bs.max() || bs.min() > f.cuts.back()){
+ int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
+ double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
+ return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));
+ }
+
+ std::vector<double> levels;//we can forget first and last cuts...
+ levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);
+ //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.
+ std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);
+
+ //-- Compose each piece of g with the relevant seg of f.
+ result.cuts.push_back(0.);
+ std::map<double,unsigned>::iterator cut=cuts_pb.begin();
+ std::map<double,unsigned>::iterator next=cut; next++;
+ while(next!=cuts_pb.end()){
+ //assert(std::abs(int((*cut).second-(*next).second))<1);
+ //TODO: find a way to recover from this error? the root finder missed some root;
+ // the levels/variations of f might be too close/fast...
+ int idx = compose_findSegIdx(cut,next,levels,g);
+ double t0=(*cut).first;
+ double t1=(*next).first;
+
+ SBasis sub_g=compose(g, Linear(t0,t1));
+ sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),
+ (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);
+ result.push(compose(f[idx],sub_g),t1);
+ cut++;
+ next++;
+ }
+ return(result);
+}
+
+//TODO: add concept check for following composition functions
+template<typename T>
+Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){
+ Piecewise<T> result;
+ for(unsigned i = 0; i < g.segs.size(); i++){
+ Piecewise<T> fgi=compose(f, g.segs[i]);
+ fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));
+ result.concat(fgi);
+ }
+ return result;
+}
+
+template <typename T>
+Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}
+template <typename T>
+Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}
+
+template<typename T>
+Piecewise<T> integral(Piecewise<T> const &a) {
+ Piecewise<T> result;
+ result.segs.resize(a.segs.size());
+ result.cuts = a.cuts;
+ typename T::output_type c = a.segs[0].at0();
+ for(unsigned i = 0; i < a.segs.size(); i++){
+ result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);
+ result.segs[i]+= c-result.segs[i].at0();
+ c = result.segs[i].at1();
+ }
+ return result;
+}
+
+template<typename T>
+Piecewise<T> derivative(Piecewise<T> const &a) {
+ Piecewise<T> result;
+ result.segs.resize(a.segs.size());
+ result.cuts = a.cuts;
+ for(unsigned i = 0; i < a.segs.size(); i++){
+ result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);
+ }
+ return result;
+}
+
+std::vector<double> roots(Piecewise<SBasis> const &f);
+
+}
+
+#endif //SEEN_GEOM_PW_SB_H
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
index 136687ae87c8fa518080f5945240e65652284521..fdc1cefe55a585d2205a63fcf925b1a74a587058 100644 (file)
-#include "poly-dk-solve.h"\r
-#include <iterator>\r
-\r
-/*** implementation of the Durand-Kerner method. seems buggy*/\r
-\r
-std::complex<double> evalu(Poly const & p, std::complex<double> x) {\r
- std::complex<double> result = 0;\r
- std::complex<double> xx = 1;\r
- \r
- for(unsigned i = 0; i < p.size(); i++) {\r
- result += p[i]*xx;\r
- xx *= x;\r
- }\r
- return result;\r
-}\r
-\r
-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {\r
- std::vector<std::complex<double> > roots;\r
- const int N = ply.degree();\r
- \r
- std::complex<double> b(0.4, 0.9);\r
- std::complex<double> p = 1;\r
- for(int i = 0; i < N; i++) {\r
- roots.push_back(p);\r
- p *= b;\r
- }\r
- assert(roots.size() == ply.degree());\r
- \r
- double error = 0;\r
- int i;\r
- for( i = 0; i < 30; i++) {\r
- error = 0;\r
- for(int r_i = 0; r_i < N; r_i++) {\r
- std::complex<double> denom = 1;\r
- std::complex<double> R = roots[r_i];\r
- for(int d_i = 0; d_i < N; d_i++) {\r
- if(r_i != d_i)\r
- denom *= R-roots[d_i];\r
- }\r
- assert(norm(denom) != 0);\r
- std::complex<double> dr = evalu(ply, R)/denom;\r
- error += norm(dr);\r
- roots[r_i] = R - dr;\r
- }\r
- /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));\r
- std::cout << std::endl;*/\r
- if(error < tol)\r
- break;\r
- }\r
- //std::cout << error << ", " << i<< std::endl;\r
- return roots;\r
-}\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+#include "poly-dk-solve.h"
+#include <iterator>
+
+/*** implementation of the Durand-Kerner method. seems buggy*/
+
+std::complex<double> evalu(Poly const & p, std::complex<double> x) {
+ std::complex<double> result = 0;
+ std::complex<double> xx = 1;
+
+ for(unsigned i = 0; i < p.size(); i++) {
+ result += p[i]*xx;
+ xx *= x;
+ }
+ return result;
+}
+
+std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
+ std::vector<std::complex<double> > roots;
+ const int N = ply.degree();
+
+ std::complex<double> b(0.4, 0.9);
+ std::complex<double> p = 1;
+ for(int i = 0; i < N; i++) {
+ roots.push_back(p);
+ p *= b;
+ }
+ assert(roots.size() == ply.degree());
+
+ double error = 0;
+ int i;
+ for( i = 0; i < 30; i++) {
+ error = 0;
+ for(int r_i = 0; r_i < N; r_i++) {
+ std::complex<double> denom = 1;
+ std::complex<double> R = roots[r_i];
+ for(int d_i = 0; d_i < N; d_i++) {
+ if(r_i != d_i)
+ denom *= R-roots[d_i];
+ }
+ assert(norm(denom) != 0);
+ std::complex<double> dr = evalu(ply, R)/denom;
+ error += norm(dr);
+ roots[r_i] = R - dr;
+ }
+ /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
+ std::cout << std::endl;*/
+ if(error < tol)
+ break;
+ }
+ //std::cout << error << ", " << i<< std::endl;
+ return roots;
+}
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
index fb3c2075b213da13d21d01a9793d675d3e9934f5..766f16edafa1e2621e9d7e24bfa0c3da9cb2f3a1 100644 (file)
-#include "poly-laguerre-solve.h"\r
-#include <iterator>\r
-\r
-typedef std::complex<double> cdouble;\r
-\r
-cdouble laguerre_internal_complex(Poly const & p,\r
- double x0,\r
- double tol,\r
- bool & quad_root) {\r
- cdouble a = 2*tol;\r
- cdouble xk = x0;\r
- double n = p.degree();\r
- quad_root = false;\r
- const unsigned shuffle_rate = 10;\r
- static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};\r
- unsigned shuffle_counter = 0;\r
- while(std::norm(a) > (tol*tol)) {\r
- //std::cout << "xk = " << xk << std::endl;\r
- cdouble b = p.back();\r
- cdouble d = 0, f = 0;\r
- double err = abs(b);\r
- double abx = abs(xk);\r
- for(int j = p.size()-2; j >= 0; j--) {\r
- f = xk*f + d;\r
- d = xk*d + b;\r
- b = xk*b + p[j];\r
- err = abs(b) + abx*err;\r
- }\r
- \r
- err *= 1e-7; // magic epsilon for convergence, should be computed from tol\r
- \r
- cdouble px = b;\r
- if(abs(b) < err)\r
- return xk;\r
- //if(std::norm(px) < tol*tol)\r
- // return xk;\r
- cdouble G = d / px;\r
- cdouble H = G*G - f / px;\r
- \r
- //std::cout << "G = " << G << "H = " << H;\r
- cdouble radicand = (n - 1)*(n*H-G*G);\r
- //assert(radicand.real() > 0);\r
- if(radicand.real() < 0)\r
- quad_root = true;\r
- //std::cout << "radicand = " << radicand << std::endl;\r
- if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation\r
- a = - sqrt(radicand);\r
- else\r
- a = sqrt(radicand);\r
- //std::cout << "a = " << a << std::endl;\r
- a = n / (a + G);\r
- //std::cout << "a = " << a << std::endl;\r
- if(shuffle_counter % shuffle_rate == 0)\r
- ;//a *= shuffle[shuffle_counter / shuffle_rate];\r
- xk -= a;\r
- shuffle_counter++;\r
- if(shuffle_counter >= 90)\r
- break;\r
- }\r
- //std::cout << "xk = " << xk << std::endl;\r
- return xk;\r
-}\r
-\r
-double laguerre_internal(Poly const & p,\r
- Poly const & pp,\r
- Poly const & ppp,\r
- double x0,\r
- double tol,\r
- bool & quad_root) {\r
- double a = 2*tol;\r
- double xk = x0;\r
- double n = p.degree();\r
- quad_root = false;\r
- while(a*a > (tol*tol)) {\r
- //std::cout << "xk = " << xk << std::endl;\r
- double px = p(xk);\r
- if(px*px < tol*tol)\r
- return xk;\r
- double G = pp(xk) / px;\r
- double H = G*G - ppp(xk) / px;\r
- \r
- //std::cout << "G = " << G << "H = " << H;\r
- double radicand = (n - 1)*(n*H-G*G);\r
- assert(radicand > 0);\r
- //std::cout << "radicand = " << radicand << std::endl;\r
- if(G < 0) // here we try to maximise the denominator avoiding cancellation\r
- a = - sqrt(radicand);\r
- else\r
- a = sqrt(radicand);\r
- //std::cout << "a = " << a << std::endl;\r
- a = n / (a + G);\r
- //std::cout << "a = " << a << std::endl;\r
- xk -= a;\r
- }\r
- //std::cout << "xk = " << xk << std::endl;\r
- return xk;\r
-}\r
-\r
-\r
-std::vector<cdouble > \r
-laguerre(Poly p, const double tol) {\r
- std::vector<cdouble > solutions;\r
- //std::cout << "p = " << p << " = ";\r
- while(p.size() > 1)\r
- {\r
- double x0 = 0;\r
- bool quad_root = false;\r
- cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);\r
- //if(abs(sol) > 1) break;\r
- Poly dvs;\r
- if(quad_root) {\r
- dvs.push_back((sol*conj(sol)).real());\r
- dvs.push_back(-(sol + conj(sol)).real());\r
- dvs.push_back(1.0);\r
- //std::cout << "(" << dvs << ")";\r
- //solutions.push_back(sol);\r
- //solutions.push_back(conj(sol));\r
- } else {\r
- //std::cout << sol << std::endl;\r
- dvs.push_back(-sol.real());\r
- dvs.push_back(1.0);\r
- solutions.push_back(sol);\r
- //std::cout << "(" << dvs << ")";\r
- }\r
- Poly r;\r
- p = divide(p, dvs, r);\r
- //std::cout << r << std::endl;\r
- }\r
- return solutions;\r
-}\r
-\r
-std::vector<double> \r
-laguerre_real_interval(Poly const & ply, \r
- const double lo, const double hi,\r
- const double tol) {\r
-}\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+#include "poly-laguerre-solve.h"
+#include <iterator>
+
+typedef std::complex<double> cdouble;
+
+cdouble laguerre_internal_complex(Poly const & p,
+ double x0,
+ double tol,
+ bool & quad_root) {
+ cdouble a = 2*tol;
+ cdouble xk = x0;
+ double n = p.degree();
+ quad_root = false;
+ const unsigned shuffle_rate = 10;
+ static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
+ unsigned shuffle_counter = 0;
+ while(std::norm(a) > (tol*tol)) {
+ //std::cout << "xk = " << xk << std::endl;
+ cdouble b = p.back();
+ cdouble d = 0, f = 0;
+ double err = abs(b);
+ double abx = abs(xk);
+ for(int j = p.size()-2; j >= 0; j--) {
+ f = xk*f + d;
+ d = xk*d + b;
+ b = xk*b + p[j];
+ err = abs(b) + abx*err;
+ }
+
+ err *= 1e-7; // magic epsilon for convergence, should be computed from tol
+
+ cdouble px = b;
+ if(abs(b) < err)
+ return xk;
+ //if(std::norm(px) < tol*tol)
+ // return xk;
+ cdouble G = d / px;
+ cdouble H = G*G - f / px;
+
+ //std::cout << "G = " << G << "H = " << H;
+ cdouble radicand = (n - 1)*(n*H-G*G);
+ //assert(radicand.real() > 0);
+ if(radicand.real() < 0)
+ quad_root = true;
+ //std::cout << "radicand = " << radicand << std::endl;
+ if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
+ a = - sqrt(radicand);
+ else
+ a = sqrt(radicand);
+ //std::cout << "a = " << a << std::endl;
+ a = n / (a + G);
+ //std::cout << "a = " << a << std::endl;
+ if(shuffle_counter % shuffle_rate == 0)
+ ;//a *= shuffle[shuffle_counter / shuffle_rate];
+ xk -= a;
+ shuffle_counter++;
+ if(shuffle_counter >= 90)
+ break;
+ }
+ //std::cout << "xk = " << xk << std::endl;
+ return xk;
+}
+
+double laguerre_internal(Poly const & p,
+ Poly const & pp,
+ Poly const & ppp,
+ double x0,
+ double tol,
+ bool & quad_root) {
+ double a = 2*tol;
+ double xk = x0;
+ double n = p.degree();
+ quad_root = false;
+ while(a*a > (tol*tol)) {
+ //std::cout << "xk = " << xk << std::endl;
+ double px = p(xk);
+ if(px*px < tol*tol)
+ return xk;
+ double G = pp(xk) / px;
+ double H = G*G - ppp(xk) / px;
+
+ //std::cout << "G = " << G << "H = " << H;
+ double radicand = (n - 1)*(n*H-G*G);
+ assert(radicand > 0);
+ //std::cout << "radicand = " << radicand << std::endl;
+ if(G < 0) // here we try to maximise the denominator avoiding cancellation
+ a = - sqrt(radicand);
+ else
+ a = sqrt(radicand);
+ //std::cout << "a = " << a << std::endl;
+ a = n / (a + G);
+ //std::cout << "a = " << a << std::endl;
+ xk -= a;
+ }
+ //std::cout << "xk = " << xk << std::endl;
+ return xk;
+}
+
+
+std::vector<cdouble >
+laguerre(Poly p, const double tol) {
+ std::vector<cdouble > solutions;
+ //std::cout << "p = " << p << " = ";
+ while(p.size() > 1)
+ {
+ double x0 = 0;
+ bool quad_root = false;
+ cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
+ //if(abs(sol) > 1) break;
+ Poly dvs;
+ if(quad_root) {
+ dvs.push_back((sol*conj(sol)).real());
+ dvs.push_back(-(sol + conj(sol)).real());
+ dvs.push_back(1.0);
+ //std::cout << "(" << dvs << ")";
+ //solutions.push_back(sol);
+ //solutions.push_back(conj(sol));
+ } else {
+ //std::cout << sol << std::endl;
+ dvs.push_back(-sol.real());
+ dvs.push_back(1.0);
+ solutions.push_back(sol);
+ //std::cout << "(" << dvs << ")";
+ }
+ Poly r;
+ p = divide(p, dvs, r);
+ //std::cout << r << std::endl;
+ }
+ return solutions;
+}
+
+std::vector<double>
+laguerre_real_interval(Poly const & ply,
+ const double lo, const double hi,
+ const double tol) {
+}
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/quadtree.h b/src/2geom/quadtree.h
index ef583c0b56fdd5a5653ff08ea8cf99ed8676fe1e..62c769c8bae5b13683e1832bf87f0ec157edc5b6 100644 (file)
--- a/src/2geom/quadtree.h
+++ b/src/2geom/quadtree.h
-#include <vector>\r
+#include <vector>
#include <cassert>
class Quad{
double by0, by1;
QuadTree() : root(0), scale(1) {}
-
+
Quad* search(double x0, double y0, double x1, double y1);
void insert(double x0, double y0, double x1, double y1, int shape);
void erase(Quad *q, int shape);
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index 5bf0d2876d14c02bd444cd637d4319c92b47b626..7157bc8856313021f405feee0f10371719ff86d2 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
-/*\r
- * sbasis.cpp - S-power basis function class + supporting classes\r
- *\r
- * Authors:\r
- * Nathan Hurst <njh@mail.csse.monash.edu.au>\r
- * Michael Sloan <mgsloan@gmail.com>\r
- *\r
- * Copyright (C) 2006-2007 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-#include <cmath>\r
-\r
-#include "sbasis.h"\r
-#include "isnan.h"\r
-\r
-namespace Geom{\r
-\r
-/*** At some point we should work on tighter bounds for the error. It is clear that the error is\r
- * bounded by the L1 norm over the tail of the series, but this is very loose, leading to far too\r
- * many cubic beziers. I've changed this to be \sum _i=tail ^\infty |hat a_i| 2^-i but I have no\r
- * evidence that this is correct.\r
- */\r
-\r
-/*\r
-double SBasis::tail_error(unsigned tail) const {\r
- double err = 0, s = 1./(1<<(2*tail)); // rough\r
- for(unsigned i = tail; i < size(); i++) {\r
- err += (fabs((*this)[i][0]) + fabs((*this)[i][1]))*s;\r
- s /= 4;\r
- }\r
- return err;\r
-}\r
+/*
+ * sbasis.cpp - S-power basis function class + supporting classes
+ *
+ * Authors:
+ * Nathan Hurst <njh@mail.csse.monash.edu.au>
+ * Michael Sloan <mgsloan@gmail.com>
+ *
+ * Copyright (C) 2006-2007 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <cmath>
+
+#include "sbasis.h"
+#include "isnan.h"
+
+namespace Geom{
+
+/*** At some point we should work on tighter bounds for the error. It is clear that the error is
+ * bounded by the L1 norm over the tail of the series, but this is very loose, leading to far too
+ * many cubic beziers. I've changed this to be \sum _i=tail ^\infty |hat a_i| 2^-i but I have no
+ * evidence that this is correct.
+ */
+
+/*
+double SBasis::tail_error(unsigned tail) const {
+ double err = 0, s = 1./(1<<(2*tail)); // rough
+ for(unsigned i = tail; i < size(); i++) {
+ err += (fabs((*this)[i][0]) + fabs((*this)[i][1]))*s;
+ s /= 4;
+ }
+ return err;
+}
*/
-\r
-double SBasis::tailError(unsigned tail) const {\r
- Interval bs = bounds_fast(*this, tail);\r
- return std::max(fabs(bs.min()),fabs(bs.max()));\r
-}\r
-\r
-bool SBasis::isFinite() const {\r
- for(unsigned i = 0; i < size(); i++) {\r
- if(!(*this)[i].isFinite())\r
- return false;\r
- }\r
- return true;\r
-}\r
-\r
-SBasis operator+(const SBasis& a, const SBasis& b) {\r
- SBasis result;\r
- const unsigned out_size = std::max(a.size(), b.size());\r
- const unsigned min_size = std::min(a.size(), b.size());\r
- result.reserve(out_size);\r
- \r
- for(unsigned i = 0; i < min_size; i++) {\r
- result.push_back(a[i] + b[i]);\r
- }\r
- for(unsigned i = min_size; i < a.size(); i++)\r
- result.push_back(a[i]);\r
- for(unsigned i = min_size; i < b.size(); i++)\r
- result.push_back(b[i]);\r
-\r
- assert(result.size() == out_size);\r
- return result;\r
-}\r
-\r
-SBasis operator-(const SBasis& a, const SBasis& b) {\r
- SBasis result;\r
- const unsigned out_size = std::max(a.size(), b.size());\r
- const unsigned min_size = std::min(a.size(), b.size());\r
- result.reserve(out_size);\r
- \r
- for(unsigned i = 0; i < min_size; i++) {\r
- result.push_back(a[i] - b[i]);\r
- }\r
- for(unsigned i = min_size; i < a.size(); i++)\r
- result.push_back(a[i]);\r
- for(unsigned i = min_size; i < b.size(); i++)\r
- result.push_back(-b[i]);\r
-\r
- assert(result.size() == out_size);\r
- return result;\r
-}\r
-\r
-SBasis& operator+=(SBasis& a, const SBasis& b) {\r
- const unsigned out_size = std::max(a.size(), b.size());\r
- const unsigned min_size = std::min(a.size(), b.size());\r
- a.reserve(out_size);\r
- \r
- for(unsigned i = 0; i < min_size; i++)\r
- a[i] += b[i];\r
- for(unsigned i = min_size; i < b.size(); i++)\r
- a.push_back(b[i]);\r
- \r
- assert(a.size() == out_size);\r
- return a;\r
-}\r
-\r
-SBasis& operator-=(SBasis& a, const SBasis& b) {\r
- const unsigned out_size = std::max(a.size(), b.size());\r
- const unsigned min_size = std::min(a.size(), b.size());\r
- a.reserve(out_size);\r
- \r
- for(unsigned i = 0; i < min_size; i++)\r
- a[i] -= b[i];\r
- for(unsigned i = min_size; i < b.size(); i++)\r
- a.push_back(-b[i]);\r
- \r
- assert(a.size() == out_size);\r
- return a;\r
-}\r
-\r
-SBasis operator*(SBasis const &a, double k) {\r
- SBasis c;\r
- c.reserve(a.size());\r
- for(unsigned i = 0; i < a.size(); i++)\r
- c.push_back(a[i] * k);\r
- return c;\r
-}\r
-\r
-SBasis& operator*=(SBasis& a, double b) {\r
- if (a.isZero()) return a;\r
- if (b == 0)\r
- a.clear();\r
- else\r
- for(unsigned i = 0; i < a.size(); i++)\r
- a[i] *= b;\r
- return a;\r
-}\r
-\r
-SBasis shift(SBasis const &a, int sh) {\r
- SBasis c = a;\r
- if(sh > 0) {\r
- c.insert(c.begin(), sh, Linear(0,0));\r
- } else {\r
- //TODO: truncate\r
- }\r
- return c;\r
-}\r
-\r
-SBasis shift(Linear const &a, int sh) {\r
- SBasis c;\r
- if(sh > 0) {\r
- c.insert(c.begin(), sh, Linear(0,0));\r
- c.push_back(a);\r
- }\r
- return c;\r
-}\r
-\r
-SBasis multiply(SBasis const &a, SBasis const &b) {\r
- // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}\r
- \r
- // shift(1, a.Tri*b.Tri)\r
- SBasis c;\r
- if(a.isZero() || b.isZero())\r
- return c;\r
- c.resize(a.size() + b.size(), Linear(0,0));\r
- c[0] = Linear(0,0);\r
- for(unsigned j = 0; j < b.size(); j++) {\r
- for(unsigned i = j; i < a.size()+j; i++) {\r
- double tri = Tri(b[j])*Tri(a[i-j]);\r
- c[i+1/*shift*/] += Linear(Hat(-tri));\r
- }\r
- }\r
- for(unsigned j = 0; j < b.size(); j++) {\r
- for(unsigned i = j; i < a.size()+j; i++) {\r
- for(unsigned dim = 0; dim < 2; dim++)\r
- c[i][dim] += b[j][dim]*a[i-j][dim];\r
- }\r
- }\r
- c.normalize();\r
- //assert(!(0 == c.back()[0] && 0 == c.back()[1]));\r
- return c;\r
-}\r
-\r
-SBasis integral(SBasis const &c) {\r
- SBasis a;\r
- a.resize(c.size() + 1, Linear(0,0));\r
- a[0] = Linear(0,0);\r
- \r
- for(unsigned k = 1; k < c.size() + 1; k++) {\r
- double ahat = -Tri(c[k-1])/(2*k);\r
- a[k] = Hat(ahat);\r
- }\r
- double aTri = 0;\r
- for(int k = c.size()-1; k >= 0; k--) {\r
- aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1);\r
- a[k][0] -= aTri/2;\r
- a[k][1] += aTri/2;\r
- }\r
- a.normalize();\r
- return a;\r
-}\r
-\r
-SBasis derivative(SBasis const &a) {\r
- SBasis c;\r
- c.resize(a.size(), Linear(0,0));\r
- \r
- for(unsigned k = 0; k < a.size(); k++) {\r
- double d = (2*k+1)*Tri(a[k]);\r
- \r
- for(unsigned dim = 0; dim < 2; dim++) {\r
- c[k][dim] = d;\r
- if(k+1 < a.size()) {\r
- if(dim)\r
- c[k][dim] = d - (k+1)*a[k+1][dim];\r
- else\r
- c[k][dim] = d + (k+1)*a[k+1][dim];\r
- }\r
- }\r
- }\r
- \r
- return c;\r
-}\r
-\r
-//TODO: convert int k to unsigned k, and remove cast\r
-SBasis sqrt(SBasis const &a, int k) {\r
- SBasis c;\r
- if(a.isZero() || k == 0)\r
- return c;\r
- c.resize(k, Linear(0,0));\r
- c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1]));\r
- SBasis r = a - multiply(c, c); // remainder\r
- \r
- for(unsigned i = 1; i <= (unsigned)k and i<r.size(); i++) {\r
- Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));\r
- SBasis cisi = shift(ci, i);\r
- r -= multiply(shift((c*2 + cisi), i), SBasis(ci));\r
- r.truncate(k+1);\r
- c += cisi;\r
- if(r.tailError(i) == 0) // if exact\r
- break;\r
- }\r
- \r
- return c;\r
-}\r
-\r
-// return a kth order approx to 1/a)\r
-SBasis reciprocal(Linear const &a, int k) {\r
- SBasis c;\r
- assert(!a.isZero());\r
- c.resize(k, Linear(0,0));\r
- double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]);\r
- double r_s0k = 1;\r
- for(unsigned i = 0; i < (unsigned)k; i++) {\r
- c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);\r
- r_s0k *= r_s0;\r
- }\r
- return c;\r
-}\r
-\r
-SBasis divide(SBasis const &a, SBasis const &b, int k) {\r
- SBasis c;\r
- assert(!a.isZero());\r
- SBasis r = a; // remainder\r
- \r
- k++;\r
- r.resize(k, Linear(0,0));\r
- c.resize(k, Linear(0,0));\r
-\r
- for(unsigned i = 0; i < (unsigned)k; i++) {\r
- Linear ci(r[i][0]/b[0][0], r[i][1]/b[0][1]); //H0\r
- c[i] += ci;\r
- r -= shift(multiply(ci,b), i);\r
- r.truncate(k+1);\r
- if(r.tailError(i) == 0) // if exact\r
- break;\r
- }\r
- \r
- return c;\r
-}\r
-\r
-// a(b)\r
-// return a0 + s(a1 + s(a2 +... where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k\r
-SBasis compose(SBasis const &a, SBasis const &b) {\r
- SBasis s = multiply((SBasis(Linear(1,1))-b), b);\r
- SBasis r;\r
- \r
- for(int i = a.size()-1; i >= 0; i--) {\r
- r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);\r
- }\r
- return r;\r
-}\r
-\r
-// a(b)\r
-// return a0 + s(a1 + s(a2 +... where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k\r
-SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {\r
- SBasis s = multiply((SBasis(Linear(1,1))-b), b);\r
- SBasis r;\r
- \r
- for(int i = a.size()-1; i >= 0; i--) {\r
- r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);\r
- }\r
- r.truncate(k);\r
- return r;\r
-}\r
-\r
-/*\r
-Inversion algorithm. The notation is certainly very misleading. The\r
-pseudocode should say:\r
-\r
-c(v) := 0\r
-r(u) := r_0(u) := u\r
-for i:=0 to k do\r
- c_i(v) := H_0(r_i(u)/(t_1)^i; u)\r
- c(v) := c(v) + c_i(v)*t^i\r
- r(u) := r(u) ? c_i(u)*(t(u))^i\r
-endfor\r
-*/\r
-\r
-//#define DEBUG_INVERSION 1 \r
-\r
-SBasis inverse(SBasis a, int k) {\r
- assert(a.size() > 0);\r
-// the function should have 'unit range'("a00 = 0 and a01 = 1") and be monotonic.\r
- double a0 = a[0][0];\r
- if(a0 != 0) {\r
- a -= a0;\r
- }\r
- double a1 = a[0][1];\r
- assert(a1 != 0); // not invertable.\r
- \r
- if(a1 != 1) {\r
- a /= a1;\r
- }\r
- SBasis c; // c(v) := 0\r
- if(a.size() >= 2 && k == 2) {\r
- c.push_back(Linear(0,1));\r
- Linear t1(1+a[1][0], 1-a[1][1]); // t_1\r
- c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]));\r
- } else if(a.size() >= 2) { // non linear\r
- SBasis r = Linear(0,1); // r(u) := r_0(u) := u\r
- Linear t1(1./(1+a[1][0]), 1./(1-a[1][1])); // 1./t_1\r
- Linear one(1,1);\r
- Linear t1i = one; // t_1^0\r
- SBasis one_minus_a = SBasis(one) - a;\r
- SBasis t = multiply(one_minus_a, a); // t(u)\r
- SBasis ti(one); // t(u)^0\r
-#ifdef DEBUG_INVERSION\r
- std::cout << "a=" << a << std::endl;\r
- std::cout << "1-a=" << one_minus_a << std::endl;\r
- std::cout << "t1=" << t1 << std::endl;\r
- //assert(t1 == t[1]);\r
-#endif\r
- \r
- c.resize(k+1, Linear(0,0));\r
- for(unsigned i = 0; i < (unsigned)k; i++) { // for i:=0 to k do\r
-#ifdef DEBUG_INVERSION\r
- std::cout << "-------" << i << ": ---------" <<std::endl;\r
- std::cout << "r=" << r << std::endl\r
- << "c=" << c << std::endl\r
- << "ti=" << ti << std::endl\r
- << std::endl;\r
-#endif\r
- if(r.size() <= i) // ensure enough space in the remainder, probably not needed\r
- r.resize(i+1, Linear(0,0));\r
- Linear ci(r[i][0]*t1i[0], r[i][1]*t1i[1]); // c_i(v) := H_0(r_i(u)/(t_1)^i; u)\r
-#ifdef DEBUG_INVERSION\r
- std::cout << "t1i=" << t1i << std::endl;\r
- std::cout << "ci=" << ci << std::endl;\r
-#endif\r
- for(int dim = 0; dim < 2; dim++) // t1^-i *= 1./t1\r
- t1i[dim] *= t1[dim];\r
- c[i] = ci; // c(v) := c(v) + c_i(v)*t^i\r
- // change from v to u parameterisation\r
- SBasis civ = one_minus_a*ci[0] + a*ci[1]; \r
- // r(u) := r(u) - c_i(u)*(t(u))^i\r
- // We can truncate this to the number of final terms, as no following terms can\r
- // contribute to the result.\r
- r -= multiply(civ,ti);\r
- r.truncate(k);\r
- if(r.tailError(i) == 0)\r
- break; // yay!\r
- ti = multiply(ti,t);\r
- }\r
-#ifdef DEBUG_INVERSION\r
- std::cout << "##########################" << std::endl;\r
-#endif\r
- } else\r
- c = Linear(0,1); // linear\r
- c -= a0; // invert the offset\r
- c /= a1; // invert the slope\r
- return c;\r
-}\r
-\r
-SBasis sin(Linear b, int k) {\r
- SBasis s = Linear(std::sin(b[0]), std::sin(b[1]));\r
- Tri tr(s[0]);\r
- double t2 = Tri(b);\r
- s.push_back(Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr));\r
- \r
- t2 *= t2;\r
- for(int i = 0; i < k; i++) {\r
- Linear bo(4*(i+1)*s[i+1][0] - 2*s[i+1][1],\r
- -2*s[i+1][0] + 4*(i+1)*s[i+1][1]);\r
- bo -= s[i]*(t2/(i+1));\r
- \r
- \r
- s.push_back(bo/double(i+2));\r
- }\r
- \r
- return s;\r
-}\r
-\r
-SBasis cos(Linear bo, int k) {\r
- return sin(Linear(bo[0] + M_PI/2,\r
- bo[1] + M_PI/2),\r
- k);\r
-}\r
-\r
-//compute fog^-1. ("zero" = double comparison threshold. *!*we might divide by "zero"*!*)\r
-//TODO: compute order according to tol?\r
-//TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!\r
-SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){\r
- SBasis result; //result\r
- SBasis r=f; //remainder\r
- SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;\r
- Pk.truncate(order);\r
- Qk.truncate(order);\r
- Pk.resize(order,Linear(0.));\r
- Qk.resize(order,Linear(0.));\r
- r.resize(order,Linear(0.));\r
-\r
- int vs= valuation(sg,zero);\r
- \r
- for (unsigned k=0; k<order; k+=vs){\r
- double p10 = Pk.at(k)[0];// we have to solve the linear system:\r
- double p01 = Pk.at(k)[1];//\r
- double q10 = Qk.at(k)[0];// p10*a + q10*b = r10\r
- double q01 = Qk.at(k)[1];// & \r
- double r10 = r.at(k)[0];// p01*a + q01*b = r01\r
- double r01 = r.at(k)[1];//\r
- double a,b;\r
- double det = p10*q01-p01*q10;\r
-\r
- //TODO: handle det~0!! \r
- if (fabs(det)<zero){\r
- det = zero;\r
- a=b=0;\r
- }else{\r
- a=( q01*r10-q10*r01)/det;\r
- b=(-p01*r10+p10*r01)/det;\r
- }\r
- result.push_back(Linear(a,b));\r
- r=r-Pk*a-Qk*b;\r
- \r
- Pk=Pk*sg;\r
- Qk=Qk*sg;\r
- Pk.truncate(order);\r
- Qk.truncate(order);\r
- r.truncate(order);\r
- }\r
- result.normalize();\r
- return result;\r
-}\r
-\r
-}\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+
+double SBasis::tailError(unsigned tail) const {
+ Interval bs = bounds_fast(*this, tail);
+ return std::max(fabs(bs.min()),fabs(bs.max()));
+}
+
+bool SBasis::isFinite() const {
+ for(unsigned i = 0; i < size(); i++) {
+ if(!(*this)[i].isFinite())
+ return false;
+ }
+ return true;
+}
+
+SBasis operator+(const SBasis& a, const SBasis& b) {
+ SBasis result;
+ const unsigned out_size = std::max(a.size(), b.size());
+ const unsigned min_size = std::min(a.size(), b.size());
+ result.reserve(out_size);
+
+ for(unsigned i = 0; i < min_size; i++) {
+ result.push_back(a[i] + b[i]);
+ }
+ for(unsigned i = min_size; i < a.size(); i++)
+ result.push_back(a[i]);
+ for(unsigned i = min_size; i < b.size(); i++)
+ result.push_back(b[i]);
+
+ assert(result.size() == out_size);
+ return result;
+}
+
+SBasis operator-(const SBasis& a, const SBasis& b) {
+ SBasis result;
+ const unsigned out_size = std::max(a.size(), b.size());
+ const unsigned min_size = std::min(a.size(), b.size());
+ result.reserve(out_size);
+
+ for(unsigned i = 0; i < min_size; i++) {
+ result.push_back(a[i] - b[i]);
+ }
+ for(unsigned i = min_size; i < a.size(); i++)
+ result.push_back(a[i]);
+ for(unsigned i = min_size; i < b.size(); i++)
+ result.push_back(-b[i]);
+
+ assert(result.size() == out_size);
+ return result;
+}
+
+SBasis& operator+=(SBasis& a, const SBasis& b) {
+ const unsigned out_size = std::max(a.size(), b.size());
+ const unsigned min_size = std::min(a.size(), b.size());
+ a.reserve(out_size);
+
+ for(unsigned i = 0; i < min_size; i++)
+ a[i] += b[i];
+ for(unsigned i = min_size; i < b.size(); i++)
+ a.push_back(b[i]);
+
+ assert(a.size() == out_size);
+ return a;
+}
+
+SBasis& operator-=(SBasis& a, const SBasis& b) {
+ const unsigned out_size = std::max(a.size(), b.size());
+ const unsigned min_size = std::min(a.size(), b.size());
+ a.reserve(out_size);
+
+ for(unsigned i = 0; i < min_size; i++)
+ a[i] -= b[i];
+ for(unsigned i = min_size; i < b.size(); i++)
+ a.push_back(-b[i]);
+
+ assert(a.size() == out_size);
+ return a;
+}
+
+SBasis operator*(SBasis const &a, double k) {
+ SBasis c;
+ c.reserve(a.size());
+ for(unsigned i = 0; i < a.size(); i++)
+ c.push_back(a[i] * k);
+ return c;
+}
+
+SBasis& operator*=(SBasis& a, double b) {
+ if (a.isZero()) return a;
+ if (b == 0)
+ a.clear();
+ else
+ for(unsigned i = 0; i < a.size(); i++)
+ a[i] *= b;
+ return a;
+}
+
+SBasis shift(SBasis const &a, int sh) {
+ SBasis c = a;
+ if(sh > 0) {
+ c.insert(c.begin(), sh, Linear(0,0));
+ } else {
+ //TODO: truncate
+ }
+ return c;
+}
+
+SBasis shift(Linear const &a, int sh) {
+ SBasis c;
+ if(sh > 0) {
+ c.insert(c.begin(), sh, Linear(0,0));
+ c.push_back(a);
+ }
+ return c;
+}
+
+SBasis multiply(SBasis const &a, SBasis const &b) {
+ // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}
+
+ // shift(1, a.Tri*b.Tri)
+ SBasis c;
+ if(a.isZero() || b.isZero())
+ return c;
+ c.resize(a.size() + b.size(), Linear(0,0));
+ c[0] = Linear(0,0);
+ for(unsigned j = 0; j < b.size(); j++) {
+ for(unsigned i = j; i < a.size()+j; i++) {
+ double tri = Tri(b[j])*Tri(a[i-j]);
+ c[i+1/*shift*/] += Linear(Hat(-tri));
+ }
+ }
+ for(unsigned j = 0; j < b.size(); j++) {
+ for(unsigned i = j; i < a.size()+j; i++) {
+ for(unsigned dim = 0; dim < 2; dim++)
+ c[i][dim] += b[j][dim]*a[i-j][dim];
+ }
+ }
+ c.normalize();
+ //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
+ return c;
+}
+
+SBasis integral(SBasis const &c) {
+ SBasis a;
+ a.resize(c.size() + 1, Linear(0,0));
+ a[0] = Linear(0,0);
+
+ for(unsigned k = 1; k < c.size() + 1; k++) {
+ double ahat = -Tri(c[k-1])/(2*k);
+ a[k] = Hat(ahat);
+ }
+ double aTri = 0;
+ for(int k = c.size()-1; k >= 0; k--) {
+ aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1);
+ a[k][0] -= aTri/2;
+ a[k][1] += aTri/2;
+ }
+ a.normalize();
+ return a;
+}
+
+SBasis derivative(SBasis const &a) {
+ SBasis c;
+ c.resize(a.size(), Linear(0,0));
+
+ for(unsigned k = 0; k < a.size(); k++) {
+ double d = (2*k+1)*Tri(a[k]);
+
+ for(unsigned dim = 0; dim < 2; dim++) {
+ c[k][dim] = d;
+ if(k+1 < a.size()) {
+ if(dim)
+ c[k][dim] = d - (k+1)*a[k+1][dim];
+ else
+ c[k][dim] = d + (k+1)*a[k+1][dim];
+ }
+ }
+ }
+
+ return c;
+}
+
+//TODO: convert int k to unsigned k, and remove cast
+SBasis sqrt(SBasis const &a, int k) {
+ SBasis c;
+ if(a.isZero() || k == 0)
+ return c;
+ c.resize(k, Linear(0,0));
+ c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1]));
+ SBasis r = a - multiply(c, c); // remainder
+
+ for(unsigned i = 1; i <= (unsigned)k and i<r.size(); i++) {
+ Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));
+ SBasis cisi = shift(ci, i);
+ r -= multiply(shift((c*2 + cisi), i), SBasis(ci));
+ r.truncate(k+1);
+ c += cisi;
+ if(r.tailError(i) == 0) // if exact
+ break;
+ }
+
+ return c;
+}
+
+// return a kth order approx to 1/a)
+SBasis reciprocal(Linear const &a, int k) {
+ SBasis c;
+ assert(!a.isZero());
+ c.resize(k, Linear(0,0));
+ double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]);
+ double r_s0k = 1;
+ for(unsigned i = 0; i < (unsigned)k; i++) {
+ c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);
+ r_s0k *= r_s0;
+ }
+ return c;
+}
+
+SBasis divide(SBasis const &a, SBasis const &b, int k) {
+ SBasis c;
+ assert(!a.isZero());
+ SBasis r = a; // remainder
+
+ k++;
+ r.resize(k, Linear(0,0));
+ c.resize(k, Linear(0,0));
+
+ for(unsigned i = 0; i < (unsigned)k; i++) {
+ Linear ci(r[i][0]/b[0][0], r[i][1]/b[0][1]); //H0
+ c[i] += ci;
+ r -= shift(multiply(ci,b), i);
+ r.truncate(k+1);
+ if(r.tailError(i) == 0) // if exact
+ break;
+ }
+
+ return c;
+}
+
+// a(b)
+// return a0 + s(a1 + s(a2 +... where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
+SBasis compose(SBasis const &a, SBasis const &b) {
+ SBasis s = multiply((SBasis(Linear(1,1))-b), b);
+ SBasis r;
+
+ for(int i = a.size()-1; i >= 0; i--) {
+ r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);
+ }
+ return r;
+}
+
+// a(b)
+// return a0 + s(a1 + s(a2 +... where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
+SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {
+ SBasis s = multiply((SBasis(Linear(1,1))-b), b);
+ SBasis r;
+
+ for(int i = a.size()-1; i >= 0; i--) {
+ r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);
+ }
+ r.truncate(k);
+ return r;
+}
+
+/*
+Inversion algorithm. The notation is certainly very misleading. The
+pseudocode should say:
+
+c(v) := 0
+r(u) := r_0(u) := u
+for i:=0 to k do
+ c_i(v) := H_0(r_i(u)/(t_1)^i; u)
+ c(v) := c(v) + c_i(v)*t^i
+ r(u) := r(u) ? c_i(u)*(t(u))^i
+endfor
+*/
+
+//#define DEBUG_INVERSION 1
+
+SBasis inverse(SBasis a, int k) {
+ assert(a.size() > 0);
+// the function should have 'unit range'("a00 = 0 and a01 = 1") and be monotonic.
+ double a0 = a[0][0];
+ if(a0 != 0) {
+ a -= a0;
+ }
+ double a1 = a[0][1];
+ assert(a1 != 0); // not invertable.
+
+ if(a1 != 1) {
+ a /= a1;
+ }
+ SBasis c; // c(v) := 0
+ if(a.size() >= 2 && k == 2) {
+ c.push_back(Linear(0,1));
+ Linear t1(1+a[1][0], 1-a[1][1]); // t_1
+ c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]));
+ } else if(a.size() >= 2) { // non linear
+ SBasis r = Linear(0,1); // r(u) := r_0(u) := u
+ Linear t1(1./(1+a[1][0]), 1./(1-a[1][1])); // 1./t_1
+ Linear one(1,1);
+ Linear t1i = one; // t_1^0
+ SBasis one_minus_a = SBasis(one) - a;
+ SBasis t = multiply(one_minus_a, a); // t(u)
+ SBasis ti(one); // t(u)^0
+#ifdef DEBUG_INVERSION
+ std::cout << "a=" << a << std::endl;
+ std::cout << "1-a=" << one_minus_a << std::endl;
+ std::cout << "t1=" << t1 << std::endl;
+ //assert(t1 == t[1]);
+#endif
+
+ c.resize(k+1, Linear(0,0));
+ for(unsigned i = 0; i < (unsigned)k; i++) { // for i:=0 to k do
+#ifdef DEBUG_INVERSION
+ std::cout << "-------" << i << ": ---------" <<std::endl;
+ std::cout << "r=" << r << std::endl
+ << "c=" << c << std::endl
+ << "ti=" << ti << std::endl
+ << std::endl;
+#endif
+ if(r.size() <= i) // ensure enough space in the remainder, probably not needed
+ r.resize(i+1, Linear(0,0));
+ Linear ci(r[i][0]*t1i[0], r[i][1]*t1i[1]); // c_i(v) := H_0(r_i(u)/(t_1)^i; u)
+#ifdef DEBUG_INVERSION
+ std::cout << "t1i=" << t1i << std::endl;
+ std::cout << "ci=" << ci << std::endl;
+#endif
+ for(int dim = 0; dim < 2; dim++) // t1^-i *= 1./t1
+ t1i[dim] *= t1[dim];
+ c[i] = ci; // c(v) := c(v) + c_i(v)*t^i
+ // change from v to u parameterisation
+ SBasis civ = one_minus_a*ci[0] + a*ci[1];
+ // r(u) := r(u) - c_i(u)*(t(u))^i
+ // We can truncate this to the number of final terms, as no following terms can
+ // contribute to the result.
+ r -= multiply(civ,ti);
+ r.truncate(k);
+ if(r.tailError(i) == 0)
+ break; // yay!
+ ti = multiply(ti,t);
+ }
+#ifdef DEBUG_INVERSION
+ std::cout << "##########################" << std::endl;
+#endif
+ } else
+ c = Linear(0,1); // linear
+ c -= a0; // invert the offset
+ c /= a1; // invert the slope
+ return c;
+}
+
+SBasis sin(Linear b, int k) {
+ SBasis s = Linear(std::sin(b[0]), std::sin(b[1]));
+ Tri tr(s[0]);
+ double t2 = Tri(b);
+ s.push_back(Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr));
+
+ t2 *= t2;
+ for(int i = 0; i < k; i++) {
+ Linear bo(4*(i+1)*s[i+1][0] - 2*s[i+1][1],
+ -2*s[i+1][0] + 4*(i+1)*s[i+1][1]);
+ bo -= s[i]*(t2/(i+1));
+
+
+ s.push_back(bo/double(i+2));
+ }
+
+ return s;
+}
+
+SBasis cos(Linear bo, int k) {
+ return sin(Linear(bo[0] + M_PI/2,
+ bo[1] + M_PI/2),
+ k);
+}
+
+//compute fog^-1. ("zero" = double comparison threshold. *!*we might divide by "zero"*!*)
+//TODO: compute order according to tol?
+//TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!
+SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){
+ SBasis result; //result
+ SBasis r=f; //remainder
+ SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;
+ Pk.truncate(order);
+ Qk.truncate(order);
+ Pk.resize(order,Linear(0.));
+ Qk.resize(order,Linear(0.));
+ r.resize(order,Linear(0.));
+
+ int vs= valuation(sg,zero);
+
+ for (unsigned k=0; k<order; k+=vs){
+ double p10 = Pk.at(k)[0];// we have to solve the linear system:
+ double p01 = Pk.at(k)[1];//
+ double q10 = Qk.at(k)[0];// p10*a + q10*b = r10
+ double q01 = Qk.at(k)[1];// &
+ double r10 = r.at(k)[0];// p01*a + q01*b = r01
+ double r01 = r.at(k)[1];//
+ double a,b;
+ double det = p10*q01-p01*q10;
+
+ //TODO: handle det~0!!
+ if (fabs(det)<zero){
+ det = zero;
+ a=b=0;
+ }else{
+ a=( q01*r10-q10*r01)/det;
+ b=(-p01*r10+p10*r01)/det;
+ }
+ result.push_back(Linear(a,b));
+ r=r-Pk*a-Qk*b;
+
+ Pk=Pk*sg;
+ Qk=Qk*sg;
+ Pk.truncate(order);
+ Qk.truncate(order);
+ r.truncate(order);
+ }
+ result.normalize();
+ return result;
+}
+
+}
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index 4329bc55925ea8eaf055672807d13121bd94f169..08674ab2fa872395a2c69da985b4e101a018e4d0 100644 (file)
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
-#include "sweep.h"\r
-\r
-#include <algorithm>\r
-\r
-namespace Geom {\r
-\r
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {\r
- std::vector<Event> events; events.reserve(rs.size()*2);\r
- std::vector<std::vector<unsigned> > pairs(rs.size());\r
- \r
- for(unsigned i = 0; i < rs.size(); i++) {\r
- events.push_back(Event(rs[i].left(), i, false));\r
- events.push_back(Event(rs[i].right(), i, true));\r
- }\r
- std::sort(events.begin(), events.end());\r
-\r
- std::vector<unsigned> open;\r
- for(unsigned i = 0; i < events.size(); i++) {\r
- unsigned ix = events[i].ix;\r
- if(events[i].closing) {\r
- std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);\r
- //if(iter != open.end())\r
- open.erase(iter);\r
- } else {\r
- for(unsigned j = 0; j < open.size(); j++) {\r
- unsigned jx = open[j];\r
- if(rs[jx][Y].intersects(rs[ix][Y])) {\r
- pairs[jx].push_back(ix);\r
- }\r
- }\r
- open.push_back(ix);\r
- }\r
- }\r
- return pairs;\r
-}\r
-\r
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {\r
- std::vector<std::vector<unsigned> > pairs(a.size());\r
- if(a.empty() || b.empty()) return pairs;\r
- std::vector<Event> events[2];\r
- events[0].reserve(a.size()*2);\r
- events[1].reserve(b.size()*2);\r
- \r
- for(unsigned n = 0; n < 2; n++) {\r
- unsigned sz = n ? b.size() : a.size();\r
- events[n].reserve(sz*2);\r
- for(unsigned i = 0; i < sz; i++) {\r
- events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));\r
- events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));\r
- }\r
- std::sort(events[n].begin(), events[n].end());\r
- }\r
-\r
- std::vector<unsigned> open[2];\r
- bool n = events[1].front() < events[0].front();\r
- for(unsigned i[] = {0,0}; i[n] < events[n].size();) {\r
- unsigned ix = events[n][i[n]].ix;\r
- bool closing = events[n][i[n]].closing;\r
- //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";\r
- if(closing) {\r
- open[n].erase(std::find(open[n].begin(), open[n].end(), ix));\r
- } else {\r
- if(n) {\r
- //n = 1\r
- //opening a B, add to all open a\r
- for(unsigned j = 0; j < open[0].size(); j++) {\r
- unsigned jx = open[0][j];\r
- if(a[jx][Y].intersects(b[ix][Y])) {\r
- pairs[jx].push_back(ix);\r
- }\r
- }\r
- } else {\r
- //n = 0\r
- //opening an A, add all open b\r
- for(unsigned j = 0; j < open[1].size(); j++) {\r
- unsigned jx = open[1][j];\r
- if(b[jx][Y].intersects(a[ix][Y])) {\r
- pairs[ix].push_back(jx);\r
- }\r
- }\r
- }\r
- open[n].push_back(ix);\r
- }\r
- i[n]++;\r
- n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;\r
- }\r
- return pairs;\r
-}\r
-\r
-//Fake cull, until the switch to the real sweep is made.\r
-std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {\r
- std::vector<std::vector<unsigned> > ret;\r
- \r
- std::vector<unsigned> all;\r
- for(unsigned j = 0; j < b; j++)\r
- all.push_back(j);\r
- \r
- for(unsigned i = 0; i < a; i++)\r
- ret.push_back(all);\r
- \r
- return ret;\r
-}\r
-\r
-}\r
+#include "sweep.h"
+
+#include <algorithm>
+
+namespace Geom {
+
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
+ std::vector<Event> events; events.reserve(rs.size()*2);
+ std::vector<std::vector<unsigned> > pairs(rs.size());
+
+ for(unsigned i = 0; i < rs.size(); i++) {
+ events.push_back(Event(rs[i].left(), i, false));
+ events.push_back(Event(rs[i].right(), i, true));
+ }
+ std::sort(events.begin(), events.end());
+
+ std::vector<unsigned> open;
+ for(unsigned i = 0; i < events.size(); i++) {
+ unsigned ix = events[i].ix;
+ if(events[i].closing) {
+ std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);
+ //if(iter != open.end())
+ open.erase(iter);
+ } else {
+ for(unsigned j = 0; j < open.size(); j++) {
+ unsigned jx = open[j];
+ if(rs[jx][Y].intersects(rs[ix][Y])) {
+ pairs[jx].push_back(ix);
+ }
+ }
+ open.push_back(ix);
+ }
+ }
+ return pairs;
+}
+
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {
+ std::vector<std::vector<unsigned> > pairs(a.size());
+ if(a.empty() || b.empty()) return pairs;
+ std::vector<Event> events[2];
+ events[0].reserve(a.size()*2);
+ events[1].reserve(b.size()*2);
+
+ for(unsigned n = 0; n < 2; n++) {
+ unsigned sz = n ? b.size() : a.size();
+ events[n].reserve(sz*2);
+ for(unsigned i = 0; i < sz; i++) {
+ events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));
+ events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));
+ }
+ std::sort(events[n].begin(), events[n].end());
+ }
+
+ std::vector<unsigned> open[2];
+ bool n = events[1].front() < events[0].front();
+ for(unsigned i[] = {0,0}; i[n] < events[n].size();) {
+ unsigned ix = events[n][i[n]].ix;
+ bool closing = events[n][i[n]].closing;
+ //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";
+ if(closing) {
+ open[n].erase(std::find(open[n].begin(), open[n].end(), ix));
+ } else {
+ if(n) {
+ //n = 1
+ //opening a B, add to all open a
+ for(unsigned j = 0; j < open[0].size(); j++) {
+ unsigned jx = open[0][j];
+ if(a[jx][Y].intersects(b[ix][Y])) {
+ pairs[jx].push_back(ix);
+ }
+ }
+ } else {
+ //n = 0
+ //opening an A, add all open b
+ for(unsigned j = 0; j < open[1].size(); j++) {
+ unsigned jx = open[1][j];
+ if(b[jx][Y].intersects(a[ix][Y])) {
+ pairs[ix].push_back(jx);
+ }
+ }
+ }
+ open[n].push_back(ix);
+ }
+ i[n]++;
+ n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;
+ }
+ return pairs;
+}
+
+//Fake cull, until the switch to the real sweep is made.
+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {
+ std::vector<std::vector<unsigned> > ret;
+
+ std::vector<unsigned> all;
+ for(unsigned j = 0; j < b; j++)
+ all.push_back(j);
+
+ for(unsigned i = 0; i < a; i++)
+ ret.push_back(all);
+
+ return ret;
+}
+
+}
index 4837317bfb903d1f7caf2cffd523e21cc89e8048..9ac75cb3cb491a29e475b31c84b8416bb48471ed 100644 (file)
-#define INKSCAPE_LPE_PATHALONGPATH_CPP\r
-\r
-/*\r
- * Copyright (C) Johan Engelen 2007 <j.b.c.engelen@utwente.nl>\r
- *\r
- * Released under GNU GPL, read the file 'COPYING' for more information\r
- */\r
-\r
-#include "live_effects/lpe-pathalongpath.h"\r
-#include "sp-shape.h"\r
-#include "sp-item.h"\r
-#include "sp-path.h"\r
-#include "display/curve.h"\r
-#include <libnr/n-art-bpath.h>\r
-#include <libnr/nr-matrix-fns.h>\r
-#include "live_effects/n-art-bpath-2geom.h"\r
-#include "svg/svg.h"\r
-#include "ui/widget/scalar.h"\r
-\r
-#include <2geom/sbasis.h>\r
-#include <2geom/sbasis-geometric.h>\r
-#include <2geom/bezier-to-sbasis.h>\r
-#include <2geom/sbasis-to-bezier.h>\r
-#include <2geom/d2.h>\r
-#include <2geom/piecewise.h>\r
-\r
-#include <algorithm>\r
-using std::vector;\r
-\r
-\r
-/* Theory in e-mail from J.F. Barraud\r
-Let B be the skeleton path, and P the pattern (the path to be deformed).\r
-\r
-P is a map t --> P(t) = ( x(t), y(t) ).\r
-B is a map t --> B(t) = ( a(t), b(t) ).\r
-\r
-The first step is to re-parametrize B by its arc length: this is the parametrization in which a point p on B is located by its distance s from start. One obtains a new map s --> U(s) = (a'(s),b'(s)), that still describes the same path B, but where the distance along B from start to\r
-U(s) is s itself.\r
-\r
-We also need a unit normal to the path. This can be obtained by computing a unit tangent vector, and rotate it by 90°. Call this normal vector N(s).\r
-\r
-The basic deformation associated to B is then given by:\r
-\r
- (x,y) --> U(x)+y*N(x)\r
-\r
-(i.e. we go for distance x along the path, and then for distance y along the normal)\r
-\r
-Of course this formula needs some minor adaptations (as is it depends on the absolute position of P for instance, so a little translation is needed\r
-first) but I think we can first forget about them.\r
-*/\r
-\r
-namespace Inkscape {\r
-namespace LivePathEffect {\r
-\r
-static const Util::EnumData<PAPCopyType> PAPCopyTypeData[PAPCT_END] = {\r
- {PAPCT_SINGLE, N_("Single"), "single"},\r
- {PAPCT_SINGLE_STRETCHED, N_("Single, stretched"), "single_stretched"},\r
- {PAPCT_REPEATED, N_("Repeated"), "repeated"},\r
- {PAPCT_REPEATED_STRETCHED, N_("Repeated, stretched"), "repeated_stretched"}\r
-};\r
-static const Util::EnumDataConverter<PAPCopyType> PAPCopyTypeConverter(PAPCopyTypeData, PAPCT_END);\r
-\r
-LPEPathAlongPath::LPEPathAlongPath(LivePathEffectObject *lpeobject) :\r
- Effect(lpeobject),\r
- bend_path(_("Bend path"), _("Path along which to bend the original path"), "bendpath", &wr, this, "M0,0 L1,0"),\r
-/* Delayed until 0.47\r
- width_path(_("Width path"), _("..."), "widthpath", &wr, this, "M0,0 L1,0"),\r
- width_path_range(_("Width path range"), _("Range of widthpath parameter"), "widthpath_range", &wr, this, 1),\r
-*/\r
- copytype(_("Path copies"), _("How many copies to place along the skeleton path"), "copytype", PAPCopyTypeConverter, &wr, this, PAPCT_SINGLE_STRETCHED),\r
- prop_scale(_("Width"), _("Width of the path"), "prop_scale", &wr, this, 1),\r
- scale_y_rel(_("Width in units of length"), _("Scale the width of the path in units of its length"), "scale_y_rel", &wr, this, false),\r
- vertical_pattern(_("Original path is vertical"), "", "vertical", &wr, this, false)\r
-{\r
- registerParameter( dynamic_cast<Parameter *>(&bend_path) );\r
-/* Delayed until 0.47\r
- registerParameter( dynamic_cast<Parameter *>(&width_path) );\r
- registerParameter( dynamic_cast<Parameter *>(&width_path_range) );\r
-*/\r
- registerParameter( dynamic_cast<Parameter *>(©type) );\r
- registerParameter( dynamic_cast<Parameter *>(&prop_scale) );\r
- registerParameter( dynamic_cast<Parameter *>(&scale_y_rel) );\r
- registerParameter( dynamic_cast<Parameter *>(&vertical_pattern) );\r
-\r
- prop_scale.param_set_digits(3);\r
- prop_scale.param_set_increments(0.01, 0.10);\r
-}\r
-\r
-LPEPathAlongPath::~LPEPathAlongPath()\r
-{\r
-\r
-}\r
-\r
-\r
-Geom::Piecewise<Geom::D2<Geom::SBasis> >\r
-LPEPathAlongPath::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > & pwd2_in)\r
-{\r
- using namespace Geom;\r
-\r
-/* Much credit should go to jfb and mgsloan of lib2geom development for the code below! */\r
-\r
- PAPCopyType type = copytype.get_value();\r
-\r
- Piecewise<D2<SBasis> > uskeleton = arc_length_parametrization(Piecewise<D2<SBasis> >(bend_path),2,.1);\r
- uskeleton = remove_short_cuts(uskeleton,.01);\r
- Rect uskeletonbounds = bounds_exact(uskeleton);\r
- uskeleton -= uskeletonbounds.midpoint();\r
-\r
- Piecewise<D2<SBasis> > n = rot90(derivative(uskeleton));\r
- n = force_continuity(remove_short_cuts(n,.1));\r
-\r
- D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pwd2_in);\r
- Piecewise<SBasis> x = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[1]) : Piecewise<SBasis>(patternd2[0]);\r
- Piecewise<SBasis> y = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[0]) : Piecewise<SBasis>(patternd2[1]);\r
- Interval pattBnds = bounds_exact(x);\r
- x -= pattBnds.min();\r
- Interval pattBndsY = bounds_exact(y);\r
- y -= pattBndsY.middle();\r
-\r
- int nbCopies = int(uskeleton.cuts.back()/pattBnds.extent());\r
- double scaling = 1;\r
-\r
- switch(type) {\r
- case PAPCT_REPEATED:\r
- break;\r
-\r
- case PAPCT_SINGLE:\r
- nbCopies = (nbCopies > 0) ? 1 : 0;\r
- break;\r
-\r
- case PAPCT_SINGLE_STRETCHED:\r
- nbCopies = 1;\r
- scaling = uskeleton.cuts.back()/pattBnds.extent();\r
- break;\r
-\r
- case PAPCT_REPEATED_STRETCHED:\r
- scaling = uskeleton.cuts.back()/(((double)nbCopies)*pattBnds.extent());\r
- break;\r
-\r
- default:\r
- return pwd2_in;\r
- };\r
-\r
- double pattWidth = pattBnds.extent() * scaling;\r
-\r
- if (scaling != 1.0) {\r
- x*=scaling;\r
- }\r
- if ( scale_y_rel.get_value() ) {\r
- y*=(scaling*prop_scale);\r
- } else {\r
- if (prop_scale != 1.0) y *= prop_scale;\r
- }\r
-\r
-/* Delayed until 0.47\r
- Piecewise<D2<SBasis> > widthpwd2 = arc_length_parametrization(Piecewise<D2<SBasis> >(width_path),2,.1);\r
- D2<Piecewise<SBasis> > widthd2pw = make_cuts_independant(widthpwd2);\r
- Piecewise<SBasis> width = (Piecewise<SBasis>(widthd2pw[Y]) - uskeletonbounds[Y].middle()) / width_path_range;\r
-*/\r
-\r
- double offs = 0;\r
- Piecewise<D2<SBasis> > output;\r
- for (int i=0; i<nbCopies; i++){\r
- output.concat(compose(uskeleton,x+offs) + y*/*compose(width,x+offs)**/compose(n,x+offs));\r
- offs+=pattWidth;\r
- }\r
-\r
- output += Point(pattBnds.middle(), pattBndsY.middle());\r
-\r
- return output;\r
-}\r
-\r
-void\r
-LPEPathAlongPath::resetDefaults(SPItem * item)\r
-{\r
- if (!SP_IS_PATH(item)) return;\r
-\r
- using namespace Geom;\r
-\r
- // set the bend path to run horizontally in the middle of the bounding box of the original path\r
- Piecewise<D2<SBasis> > pwd2;\r
- std::vector<Path> temppath = SVGD_to_2GeomPath( SP_OBJECT_REPR(item)->attribute("inkscape:original-d"));\r
- for (unsigned int i=0; i < temppath.size(); i++) {\r
- pwd2.concat( temppath[i].toPwSb() );\r
- }\r
-\r
- D2<Piecewise<SBasis> > d2pw = make_cuts_independant(pwd2);\r
- Interval bndsX = bounds_exact(d2pw[0]);\r
- Interval bndsY = bounds_exact(d2pw[1]);\r
- Point start(bndsX.min(), (bndsY.max()+bndsY.min())/2);\r
- Point end(bndsX.max(), (bndsY.max()+bndsY.min())/2);\r
-\r
- Geom::Path path;\r
- path.start( start );\r
- path.appendNew<Geom::LineSegment>( end );\r
- bend_path.param_set_and_write_new_value( path.toPwSb() );\r
-\r
-/* Delayed until 0.47\r
- Point startw(bndsX.min(), bndsY.max());\r
- Point endw(bndsX.max(), bndsY.max());\r
- Geom::Path pathw;\r
- pathw.start( startw );\r
- pathw.appendNew<Geom::LineSegment>( endw );\r
- width_path.param_set_and_write_new_value( pathw.toPwSb() );\r
- width_path_range.param_set_value(startw[Y]-start[Y]);\r
-*/\r
-}\r
-\r
-} // namespace LivePathEffect\r
-} /* namespace Inkscape */\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :\r
+#define INKSCAPE_LPE_PATHALONGPATH_CPP
+
+/*
+ * Copyright (C) Johan Engelen 2007 <j.b.c.engelen@utwente.nl>
+ *
+ * Released under GNU GPL, read the file 'COPYING' for more information
+ */
+
+#include "live_effects/lpe-pathalongpath.h"
+#include "sp-shape.h"
+#include "sp-item.h"
+#include "sp-path.h"
+#include "display/curve.h"
+#include <libnr/n-art-bpath.h>
+#include <libnr/nr-matrix-fns.h>
+#include "live_effects/n-art-bpath-2geom.h"
+#include "svg/svg.h"
+#include "ui/widget/scalar.h"
+
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-geometric.h>
+#include <2geom/bezier-to-sbasis.h>
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/d2.h>
+#include <2geom/piecewise.h>
+
+#include <algorithm>
+using std::vector;
+
+
+/* Theory in e-mail from J.F. Barraud
+Let B be the skeleton path, and P the pattern (the path to be deformed).
+
+P is a map t --> P(t) = ( x(t), y(t) ).
+B is a map t --> B(t) = ( a(t), b(t) ).
+
+The first step is to re-parametrize B by its arc length: this is the parametrization in which a point p on B is located by its distance s from start. One obtains a new map s --> U(s) = (a'(s),b'(s)), that still describes the same path B, but where the distance along B from start to
+U(s) is s itself.
+
+We also need a unit normal to the path. This can be obtained by computing a unit tangent vector, and rotate it by 90°. Call this normal vector N(s).
+
+The basic deformation associated to B is then given by:
+
+ (x,y) --> U(x)+y*N(x)
+
+(i.e. we go for distance x along the path, and then for distance y along the normal)
+
+Of course this formula needs some minor adaptations (as is it depends on the absolute position of P for instance, so a little translation is needed
+first) but I think we can first forget about them.
+*/
+
+namespace Inkscape {
+namespace LivePathEffect {
+
+static const Util::EnumData<PAPCopyType> PAPCopyTypeData[PAPCT_END] = {
+ {PAPCT_SINGLE, N_("Single"), "single"},
+ {PAPCT_SINGLE_STRETCHED, N_("Single, stretched"), "single_stretched"},
+ {PAPCT_REPEATED, N_("Repeated"), "repeated"},
+ {PAPCT_REPEATED_STRETCHED, N_("Repeated, stretched"), "repeated_stretched"}
+};
+static const Util::EnumDataConverter<PAPCopyType> PAPCopyTypeConverter(PAPCopyTypeData, PAPCT_END);
+
+LPEPathAlongPath::LPEPathAlongPath(LivePathEffectObject *lpeobject) :
+ Effect(lpeobject),
+ bend_path(_("Bend path"), _("Path along which to bend the original path"), "bendpath", &wr, this, "M0,0 L1,0"),
+/* Delayed until 0.47
+ width_path(_("Width path"), _("..."), "widthpath", &wr, this, "M0,0 L1,0"),
+ width_path_range(_("Width path range"), _("Range of widthpath parameter"), "widthpath_range", &wr, this, 1),
+*/
+ copytype(_("Path copies"), _("How many copies to place along the skeleton path"), "copytype", PAPCopyTypeConverter, &wr, this, PAPCT_SINGLE_STRETCHED),
+ prop_scale(_("Width"), _("Width of the path"), "prop_scale", &wr, this, 1),
+ scale_y_rel(_("Width in units of length"), _("Scale the width of the path in units of its length"), "scale_y_rel", &wr, this, false),
+ vertical_pattern(_("Original path is vertical"), "", "vertical", &wr, this, false)
+{
+ registerParameter( dynamic_cast<Parameter *>(&bend_path) );
+/* Delayed until 0.47
+ registerParameter( dynamic_cast<Parameter *>(&width_path) );
+ registerParameter( dynamic_cast<Parameter *>(&width_path_range) );
+*/
+ registerParameter( dynamic_cast<Parameter *>(©type) );
+ registerParameter( dynamic_cast<Parameter *>(&prop_scale) );
+ registerParameter( dynamic_cast<Parameter *>(&scale_y_rel) );
+ registerParameter( dynamic_cast<Parameter *>(&vertical_pattern) );
+
+ prop_scale.param_set_digits(3);
+ prop_scale.param_set_increments(0.01, 0.10);
+}
+
+LPEPathAlongPath::~LPEPathAlongPath()
+{
+
+}
+
+
+Geom::Piecewise<Geom::D2<Geom::SBasis> >
+LPEPathAlongPath::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > & pwd2_in)
+{
+ using namespace Geom;
+
+/* Much credit should go to jfb and mgsloan of lib2geom development for the code below! */
+
+ PAPCopyType type = copytype.get_value();
+
+ Piecewise<D2<SBasis> > uskeleton = arc_length_parametrization(Piecewise<D2<SBasis> >(bend_path),2,.1);
+ uskeleton = remove_short_cuts(uskeleton,.01);
+ Rect uskeletonbounds = bounds_exact(uskeleton);
+ uskeleton -= uskeletonbounds.midpoint();
+
+ Piecewise<D2<SBasis> > n = rot90(derivative(uskeleton));
+ n = force_continuity(remove_short_cuts(n,.1));
+
+ D2<Piecewise<SBasis> > patternd2 = make_cuts_independant(pwd2_in);
+ Piecewise<SBasis> x = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[1]) : Piecewise<SBasis>(patternd2[0]);
+ Piecewise<SBasis> y = vertical_pattern.get_value() ? Piecewise<SBasis>(patternd2[0]) : Piecewise<SBasis>(patternd2[1]);
+ Interval pattBnds = bounds_exact(x);
+ x -= pattBnds.min();
+ Interval pattBndsY = bounds_exact(y);
+ y -= pattBndsY.middle();
+
+ int nbCopies = int(uskeleton.cuts.back()/pattBnds.extent());
+ double scaling = 1;
+
+ switch(type) {
+ case PAPCT_REPEATED:
+ break;
+
+ case PAPCT_SINGLE:
+ nbCopies = (nbCopies > 0) ? 1 : 0;
+ break;
+
+ case PAPCT_SINGLE_STRETCHED:
+ nbCopies = 1;
+ scaling = uskeleton.cuts.back()/pattBnds.extent();
+ break;
+
+ case PAPCT_REPEATED_STRETCHED:
+ scaling = uskeleton.cuts.back()/(((double)nbCopies)*pattBnds.extent());
+ break;
+
+ default:
+ return pwd2_in;
+ };
+
+ double pattWidth = pattBnds.extent() * scaling;
+
+ if (scaling != 1.0) {
+ x*=scaling;
+ }
+ if ( scale_y_rel.get_value() ) {
+ y*=(scaling*prop_scale);
+ } else {
+ if (prop_scale != 1.0) y *= prop_scale;
+ }
+
+/* Delayed until 0.47
+ Piecewise<D2<SBasis> > widthpwd2 = arc_length_parametrization(Piecewise<D2<SBasis> >(width_path),2,.1);
+ D2<Piecewise<SBasis> > widthd2pw = make_cuts_independant(widthpwd2);
+ Piecewise<SBasis> width = (Piecewise<SBasis>(widthd2pw[Y]) - uskeletonbounds[Y].middle()) / width_path_range;
+*/
+
+ double offs = 0;
+ Piecewise<D2<SBasis> > output;
+ for (int i=0; i<nbCopies; i++){
+ output.concat(compose(uskeleton,x+offs) + y*/*compose(width,x+offs)**/compose(n,x+offs));
+ offs+=pattWidth;
+ }
+
+ output += Point(pattBnds.middle(), pattBndsY.middle());
+
+ return output;
+}
+
+void
+LPEPathAlongPath::resetDefaults(SPItem * item)
+{
+ if (!SP_IS_PATH(item)) return;
+
+ using namespace Geom;
+
+ // set the bend path to run horizontally in the middle of the bounding box of the original path
+ Piecewise<D2<SBasis> > pwd2;
+ std::vector<Path> temppath = SVGD_to_2GeomPath( SP_OBJECT_REPR(item)->attribute("inkscape:original-d"));
+ for (unsigned int i=0; i < temppath.size(); i++) {
+ pwd2.concat( temppath[i].toPwSb() );
+ }
+
+ D2<Piecewise<SBasis> > d2pw = make_cuts_independant(pwd2);
+ Interval bndsX = bounds_exact(d2pw[0]);
+ Interval bndsY = bounds_exact(d2pw[1]);
+ Point start(bndsX.min(), (bndsY.max()+bndsY.min())/2);
+ Point end(bndsX.max(), (bndsY.max()+bndsY.min())/2);
+
+ Geom::Path path;
+ path.start( start );
+ path.appendNew<Geom::LineSegment>( end );
+ bend_path.param_set_and_write_new_value( path.toPwSb() );
+
+/* Delayed until 0.47
+ Point startw(bndsX.min(), bndsY.max());
+ Point endw(bndsX.max(), bndsY.max());
+ Geom::Path pathw;
+ pathw.start( startw );
+ pathw.appendNew<Geom::LineSegment>( endw );
+ width_path.param_set_and_write_new_value( pathw.toPwSb() );
+ width_path_range.param_set_value(startw[Y]-start[Y]);
+*/
+}
+
+} // namespace LivePathEffect
+} /* namespace Inkscape */
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
index be064dd044b3d79d3783b646a8572b4a1da360a4..b556efd88fd30f612d33ca351b40f0de4aa144b6 100644 (file)
-#ifndef INKSCAPE_LPE_PATHALONGPATH_H\r
-#define INKSCAPE_LPE_PATHALONGPATH_H\r
-\r
-/*\r
- * Inkscape::LPEPathAlongPath\r
- *\r
-* Copyright (C) Johan Engelen 2007 <j.b.c.engelen@utwente.nl>\r
- *\r
- * Released under GNU GPL, read the file 'COPYING' for more information\r
- */\r
-\r
-#include "live_effects/effect.h"\r
-#include "live_effects/parameter/path.h"\r
-#include "live_effects/parameter/enum.h"\r
-#include "live_effects/parameter/bool.h"\r
-\r
-namespace Inkscape {\r
-namespace LivePathEffect {\r
-\r
-enum PAPCopyType {\r
- PAPCT_SINGLE = 0,\r
- PAPCT_SINGLE_STRETCHED,\r
- PAPCT_REPEATED,\r
- PAPCT_REPEATED_STRETCHED,\r
- PAPCT_END // This must be last\r
-};\r
-\r
-class LPEPathAlongPath : public Effect {\r
-public:\r
- LPEPathAlongPath(LivePathEffectObject *lpeobject);\r
- virtual ~LPEPathAlongPath();\r
-\r
- virtual Geom::Piecewise<Geom::D2<Geom::SBasis> > doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > & pwd2_in);\r
-\r
- virtual void resetDefaults(SPItem * item);\r
-\r
-private:\r
- PathParam bend_path;\r
- //PathParam width_path;\r
- //ScalarParam width_path_range;\r
- EnumParam<PAPCopyType> copytype;\r
- ScalarParam prop_scale;\r
- BoolParam scale_y_rel;\r
- BoolParam vertical_pattern;\r
-\r
- void on_pattern_pasted();\r
-\r
- LPEPathAlongPath(const LPEPathAlongPath&);\r
- LPEPathAlongPath& operator=(const LPEPathAlongPath&);\r
-};\r
-\r
-}; //namespace LivePathEffect\r
-}; //namespace Inkscape\r
-\r
-#endif\r
+#ifndef INKSCAPE_LPE_PATHALONGPATH_H
+#define INKSCAPE_LPE_PATHALONGPATH_H
+
+/*
+ * Inkscape::LPEPathAlongPath
+ *
+* Copyright (C) Johan Engelen 2007 <j.b.c.engelen@utwente.nl>
+ *
+ * Released under GNU GPL, read the file 'COPYING' for more information
+ */
+
+#include "live_effects/effect.h"
+#include "live_effects/parameter/path.h"
+#include "live_effects/parameter/enum.h"
+#include "live_effects/parameter/bool.h"
+
+namespace Inkscape {
+namespace LivePathEffect {
+
+enum PAPCopyType {
+ PAPCT_SINGLE = 0,
+ PAPCT_SINGLE_STRETCHED,
+ PAPCT_REPEATED,
+ PAPCT_REPEATED_STRETCHED,
+ PAPCT_END // This must be last
+};
+
+class LPEPathAlongPath : public Effect {
+public:
+ LPEPathAlongPath(LivePathEffectObject *lpeobject);
+ virtual ~LPEPathAlongPath();
+
+ virtual Geom::Piecewise<Geom::D2<Geom::SBasis> > doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > & pwd2_in);
+
+ virtual void resetDefaults(SPItem * item);
+
+private:
+ PathParam bend_path;
+ //PathParam width_path;
+ //ScalarParam width_path_range;
+ EnumParam<PAPCopyType> copytype;
+ ScalarParam prop_scale;
+ BoolParam scale_y_rel;
+ BoolParam vertical_pattern;
+
+ void on_pattern_pasted();
+
+ LPEPathAlongPath(const LPEPathAlongPath&);
+ LPEPathAlongPath& operator=(const LPEPathAlongPath&);
+};
+
+}; //namespace LivePathEffect
+}; //namespace Inkscape
+
+#endif