From 182486105ac34c36a2c74e775bcb44d5cbbaffac Mon Sep 17 00:00:00 2001 From: joncruz Date: Mon, 17 Dec 2007 06:59:40 +0000 Subject: [PATCH] CRLF fix --- src/2geom/ord.h | 30 +- src/2geom/piecewise.h | 1380 ++++++++++++------------ src/2geom/poly-dk-solve.cpp | 128 +-- src/2geom/poly-laguerre-solve.cpp | 294 ++--- src/2geom/quadtree.h | 4 +- src/2geom/sbasis.cpp | 978 ++++++++--------- src/2geom/sweep.cpp | 208 ++-- src/live_effects/lpe-pathalongpath.cpp | 442 ++++---- src/live_effects/lpe-pathalongpath.h | 110 +- 9 files changed, 1787 insertions(+), 1787 deletions(-) diff --git a/src/2geom/ord.h b/src/2geom/ord.h index d0e348aec..735de6a3a 100644 --- a/src/2geom/ord.h +++ b/src/2geom/ord.h @@ -4,11 +4,11 @@ namespace { -enum Cmp { - LESS_THAN=-1, - GREATER_THAN=1, - EQUAL_TO=0 -}; +enum Cmp { + LESS_THAN=-1, + GREATER_THAN=1, + EQUAL_TO=0 +}; inline Cmp operator-(Cmp x) { switch(x) { @@ -20,16 +20,16 @@ inline Cmp operator-(Cmp x) { return EQUAL_TO; } } - -template -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; - } + +template +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 c70ecd42c..171599768 100644 --- a/src/2geom/piecewise.h +++ b/src/2geom/piecewise.h @@ -1,690 +1,690 @@ -/* - * piecewise.h - Piecewise function class - * - * Copyright 2007 Michael Sloan - * - * 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 -#include - -#include "concepts.h" -#include "isnan.h" -#include - -namespace Geom { - -template -class Piecewise { - BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept); - - public: - std::vector cuts; - std::vector 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 operator()(SBasis f); - Piecewise operator()(Piecewisef); - - 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 &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 &other) { - boost::function_requires >(); - 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 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 -inline typename FragmentConcept::BoundsType bounds_fast(const Piecewise &f) { - boost::function_requires >(); - - if(f.empty()) return typename FragmentConcept::BoundsType(); - typename FragmentConcept::BoundsType ret(bounds_fast(f[0])); - for(unsigned i = 1; i < f.size(); i++) - ret.unionWith(bounds_fast(f[i])); - return ret; -} - -template -inline typename FragmentConcept::BoundsType bounds_exact(const Piecewise &f) { - boost::function_requires >(); - - if(f.empty()) return typename FragmentConcept::BoundsType(); - typename FragmentConcept::BoundsType ret(bounds_exact(f[0])); - for(unsigned i = 1; i < f.size(); i++) - ret.unionWith(bounds_exact(f[i])); - return ret; -} - -template -inline typename FragmentConcept::BoundsType bounds_local(const Piecewise &f, const Interval &m) { - boost::function_requires >(); - - if(f.empty()) return typename FragmentConcept::BoundsType(); - if(m.isEmpty()) return typename FragmentConcept::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::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, given the piece's index and a to/from time. -template -T elem_portion(const Piecewise &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 partition(const Piecewise &pw, std::vector const &c); - * Further subdivides the Piecewise such that there is a cut at every value in c. - * Precondition: c sorted lower to higher. - * - * //Given Piecewise a and b: - * Piecewise ac = a.partition(b.cuts); - * Piecewise bc = b.partition(a.cuts); - * //ac.cuts should be equivalent to bc.cuts - */ -template -Piecewise partition(const Piecewise &pw, std::vector const &c) { - assert(pw.invariants()); - if(c.empty()) return Piecewise(pw); - - Piecewise ret = Piecewise(); - 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, 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 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, 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 portion(const Piecewise &pw, double from, double to); - * Returns a Piecewise with a defined domain of [min(from, to), max(from, to)]. - */ -template -Piecewise portion(const Piecewise &pw, double from, double to) { - if(pw.empty() || from == to) return Piecewise(); - - Piecewise 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 -Piecewise remove_short_cuts(Piecewise const &f, double tol) { - if(f.empty()) return f; - Piecewise ret; - ret.push_cut(f.cuts[0]); - for(unsigned i=0; i= tol || i==f.size()-1) { - ret.push(f[i], f.cuts[i+1]); - } - } - return ret; -} - -template -Piecewise remove_short_cuts_extending(Piecewise const &f, double tol) { - if(f.empty()) return f; - Piecewise ret; - ret.push_cut(f.cuts[0]); - double last = f.cuts[0]; // last cut included - for(unsigned i=0; 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 -std::vector roots(const Piecewise &pw) { - std::vector ret; - for(unsigned i = 0; i < pw.size(); i++) { - std::vector 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 -Piecewise operator+(Piecewise const &a, typename T::output_type b) { - boost::function_requires >(); -//TODO:empty - Piecewise ret = Piecewise(); - ret.cuts = a.cuts; - for(unsigned i = 0; i < a.size();i++) - ret.push_seg(a[i] + b); - return ret; -} -template -Piecewise operator-(Piecewise const &a, typename T::output_type b) { - boost::function_requires >(); -//TODO: empty - Piecewise ret = Piecewise(); - ret.cuts = a.cuts; - for(unsigned i = 0; i < a.size();i++) - ret.push_seg(a[i] - b); - return ret; -} -template -Piecewise operator+=(Piecewise& a, typename T::output_type b) { - boost::function_requires >(); - - 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 -Piecewise operator-=(Piecewise& a, typename T::output_type b) { - boost::function_requires >(); - - 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 -Piecewise operator-(Piecewise const &a) { - boost::function_requires >(); - - Piecewise ret; - ret.cuts = a.cuts; - for(unsigned i = 0; i < a.size();i++) - ret.push_seg(- a[i]); - return ret; -} -template -Piecewise operator*(Piecewise const &a, double b) { - boost::function_requires >(); - - if(a.empty()) return Piecewise(); - - Piecewise ret; - ret.cuts = a.cuts; - for(unsigned i = 0; i < a.size();i++) - ret.push_seg(a[i] * b); - return ret; -} -template -Piecewise operator/(Piecewise const &a, double b) { - boost::function_requires >(); - - //FIXME: b == 0? - if(a.empty()) return Piecewise(); - - Piecewise ret; - ret.cuts = a.cuts; - for(unsigned i = 0; i < a.size();i++) - ret.push_seg(a[i] / b); - return ret; -} -template -Piecewise operator*=(Piecewise& a, double b) { - boost::function_requires >(); - - if(a.empty()) return Piecewise(); - - for(unsigned i = 0; i < a.size();i++) - a[i] *= b; - return a; -} -template -Piecewise operator/=(Piecewise& a, double b) { - boost::function_requires >(); - - //FIXME: b == 0? - if(a.empty()) return Piecewise(); - - for(unsigned i = 0; i < a.size();i++) - a[i] /= b; - return a; -} - -//IMPL: AddableConcept -template -Piecewise operator+(Piecewise const &a, Piecewise const &b) { - boost::function_requires >(); - - Piecewise pa = partition(a, b.cuts), pb = partition(b, a.cuts); - Piecewise ret = Piecewise(); - 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 -Piecewise operator-(Piecewise const &a, Piecewise const &b) { - boost::function_requires >(); - - Piecewise pa = partition(a, b.cuts), pb = partition(b, a.cuts); - Piecewise ret = Piecewise(); - 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 -inline Piecewise operator+=(Piecewise &a, Piecewise const &b) { - a = a+b; - return a; -} -template -inline Piecewise operator-=(Piecewise &a, Piecewise const &b) { - a = a-b; - return a; -} - -template -Piecewise operator*(Piecewise const &a, Piecewise const &b) { - //function_requires >(); - //function_requires >(); - - Piecewise pa = partition(a, b.cuts); - Piecewise pb = partition(b, a.cuts); - Piecewise ret = Piecewise(); - 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 -inline Piecewise operator*=(Piecewise &a, Piecewise const &b) { - a = a * b; - return a; -} - -Piecewise divide(Piecewise const &a, Piecewise 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 -divide(Piecewise const &a, Piecewise const &b, double tol, unsigned k, double zero=1.e-3); -Piecewise -divide(SBasis const &a, Piecewise const &b, double tol, unsigned k, double zero=1.e-3); -Piecewise -divide(Piecewise const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3); -Piecewise -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 compose_pullback(std::vector const &cuts, SBasis const &g); -int compose_findSegIdx(std::map::iterator const &cut, - std::map::iterator const &next, - std::vector const &levels, - SBasis const &g); - -//TODO: add concept check -template -Piecewise compose(Piecewise const &f, SBasis const &g){ - Piecewise result; - if (f.empty()) return result; - if (g.isZero()) return Piecewise(f(0)); - if (f.size()==1){ - double t0 = f.cuts[0], width = f.cuts[1] - t0; - return (Piecewise) 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) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g)); - } - - std::vector 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 > instead of a map. - std::map cuts_pb = compose_pullback(levels,g); - - //-- Compose each piece of g with the relevant seg of f. - result.cuts.push_back(0.); - std::map::iterator cut=cuts_pb.begin(); - std::map::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 -Piecewise compose(Piecewise const &f, Piecewise const &g){ - Piecewise result; - for(unsigned i = 0; i < g.segs.size(); i++){ - Piecewise fgi=compose(f, g.segs[i]); - fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1])); - result.concat(fgi); - } - return result; -} - -template -Piecewise Piecewise::operator()(SBasis f){return compose((*this),f);} -template -Piecewise Piecewise::operator()(Piecewisef){return compose((*this),f);} - -template -Piecewise integral(Piecewise const &a) { - Piecewise 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 -Piecewise derivative(Piecewise const &a) { - Piecewise 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 roots(Piecewise 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 : +/* + * piecewise.h - Piecewise function class + * + * Copyright 2007 Michael Sloan + * + * 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 +#include + +#include "concepts.h" +#include "isnan.h" +#include + +namespace Geom { + +template +class Piecewise { + BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept); + + public: + std::vector cuts; + std::vector 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 operator()(SBasis f); + Piecewise operator()(Piecewisef); + + 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 &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 &other) { + boost::function_requires >(); + 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 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 +inline typename FragmentConcept::BoundsType bounds_fast(const Piecewise &f) { + boost::function_requires >(); + + if(f.empty()) return typename FragmentConcept::BoundsType(); + typename FragmentConcept::BoundsType ret(bounds_fast(f[0])); + for(unsigned i = 1; i < f.size(); i++) + ret.unionWith(bounds_fast(f[i])); + return ret; +} + +template +inline typename FragmentConcept::BoundsType bounds_exact(const Piecewise &f) { + boost::function_requires >(); + + if(f.empty()) return typename FragmentConcept::BoundsType(); + typename FragmentConcept::BoundsType ret(bounds_exact(f[0])); + for(unsigned i = 1; i < f.size(); i++) + ret.unionWith(bounds_exact(f[i])); + return ret; +} + +template +inline typename FragmentConcept::BoundsType bounds_local(const Piecewise &f, const Interval &m) { + boost::function_requires >(); + + if(f.empty()) return typename FragmentConcept::BoundsType(); + if(m.isEmpty()) return typename FragmentConcept::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::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, given the piece's index and a to/from time. +template +T elem_portion(const Piecewise &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 partition(const Piecewise &pw, std::vector const &c); + * Further subdivides the Piecewise such that there is a cut at every value in c. + * Precondition: c sorted lower to higher. + * + * //Given Piecewise a and b: + * Piecewise ac = a.partition(b.cuts); + * Piecewise bc = b.partition(a.cuts); + * //ac.cuts should be equivalent to bc.cuts + */ +template +Piecewise partition(const Piecewise &pw, std::vector const &c) { + assert(pw.invariants()); + if(c.empty()) return Piecewise(pw); + + Piecewise ret = Piecewise(); + 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, 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 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, 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 portion(const Piecewise &pw, double from, double to); + * Returns a Piecewise with a defined domain of [min(from, to), max(from, to)]. + */ +template +Piecewise portion(const Piecewise &pw, double from, double to) { + if(pw.empty() || from == to) return Piecewise(); + + Piecewise 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 +Piecewise remove_short_cuts(Piecewise const &f, double tol) { + if(f.empty()) return f; + Piecewise ret; + ret.push_cut(f.cuts[0]); + for(unsigned i=0; i= tol || i==f.size()-1) { + ret.push(f[i], f.cuts[i+1]); + } + } + return ret; +} + +template +Piecewise remove_short_cuts_extending(Piecewise const &f, double tol) { + if(f.empty()) return f; + Piecewise ret; + ret.push_cut(f.cuts[0]); + double last = f.cuts[0]; // last cut included + for(unsigned i=0; 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 +std::vector roots(const Piecewise &pw) { + std::vector ret; + for(unsigned i = 0; i < pw.size(); i++) { + std::vector 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 +Piecewise operator+(Piecewise const &a, typename T::output_type b) { + boost::function_requires >(); +//TODO:empty + Piecewise ret = Piecewise(); + ret.cuts = a.cuts; + for(unsigned i = 0; i < a.size();i++) + ret.push_seg(a[i] + b); + return ret; +} +template +Piecewise operator-(Piecewise const &a, typename T::output_type b) { + boost::function_requires >(); +//TODO: empty + Piecewise ret = Piecewise(); + ret.cuts = a.cuts; + for(unsigned i = 0; i < a.size();i++) + ret.push_seg(a[i] - b); + return ret; +} +template +Piecewise operator+=(Piecewise& a, typename T::output_type b) { + boost::function_requires >(); + + 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 +Piecewise operator-=(Piecewise& a, typename T::output_type b) { + boost::function_requires >(); + + 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 +Piecewise operator-(Piecewise const &a) { + boost::function_requires >(); + + Piecewise ret; + ret.cuts = a.cuts; + for(unsigned i = 0; i < a.size();i++) + ret.push_seg(- a[i]); + return ret; +} +template +Piecewise operator*(Piecewise const &a, double b) { + boost::function_requires >(); + + if(a.empty()) return Piecewise(); + + Piecewise ret; + ret.cuts = a.cuts; + for(unsigned i = 0; i < a.size();i++) + ret.push_seg(a[i] * b); + return ret; +} +template +Piecewise operator/(Piecewise const &a, double b) { + boost::function_requires >(); + + //FIXME: b == 0? + if(a.empty()) return Piecewise(); + + Piecewise ret; + ret.cuts = a.cuts; + for(unsigned i = 0; i < a.size();i++) + ret.push_seg(a[i] / b); + return ret; +} +template +Piecewise operator*=(Piecewise& a, double b) { + boost::function_requires >(); + + if(a.empty()) return Piecewise(); + + for(unsigned i = 0; i < a.size();i++) + a[i] *= b; + return a; +} +template +Piecewise operator/=(Piecewise& a, double b) { + boost::function_requires >(); + + //FIXME: b == 0? + if(a.empty()) return Piecewise(); + + for(unsigned i = 0; i < a.size();i++) + a[i] /= b; + return a; +} + +//IMPL: AddableConcept +template +Piecewise operator+(Piecewise const &a, Piecewise const &b) { + boost::function_requires >(); + + Piecewise pa = partition(a, b.cuts), pb = partition(b, a.cuts); + Piecewise ret = Piecewise(); + 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 +Piecewise operator-(Piecewise const &a, Piecewise const &b) { + boost::function_requires >(); + + Piecewise pa = partition(a, b.cuts), pb = partition(b, a.cuts); + Piecewise ret = Piecewise(); + 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 +inline Piecewise operator+=(Piecewise &a, Piecewise const &b) { + a = a+b; + return a; +} +template +inline Piecewise operator-=(Piecewise &a, Piecewise const &b) { + a = a-b; + return a; +} + +template +Piecewise operator*(Piecewise const &a, Piecewise const &b) { + //function_requires >(); + //function_requires >(); + + Piecewise pa = partition(a, b.cuts); + Piecewise pb = partition(b, a.cuts); + Piecewise ret = Piecewise(); + 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 +inline Piecewise operator*=(Piecewise &a, Piecewise const &b) { + a = a * b; + return a; +} + +Piecewise divide(Piecewise const &a, Piecewise 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 +divide(Piecewise const &a, Piecewise const &b, double tol, unsigned k, double zero=1.e-3); +Piecewise +divide(SBasis const &a, Piecewise const &b, double tol, unsigned k, double zero=1.e-3); +Piecewise +divide(Piecewise const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3); +Piecewise +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 compose_pullback(std::vector const &cuts, SBasis const &g); +int compose_findSegIdx(std::map::iterator const &cut, + std::map::iterator const &next, + std::vector const &levels, + SBasis const &g); + +//TODO: add concept check +template +Piecewise compose(Piecewise const &f, SBasis const &g){ + Piecewise result; + if (f.empty()) return result; + if (g.isZero()) return Piecewise(f(0)); + if (f.size()==1){ + double t0 = f.cuts[0], width = f.cuts[1] - t0; + return (Piecewise) 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) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g)); + } + + std::vector 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 > instead of a map. + std::map cuts_pb = compose_pullback(levels,g); + + //-- Compose each piece of g with the relevant seg of f. + result.cuts.push_back(0.); + std::map::iterator cut=cuts_pb.begin(); + std::map::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 +Piecewise compose(Piecewise const &f, Piecewise const &g){ + Piecewise result; + for(unsigned i = 0; i < g.segs.size(); i++){ + Piecewise fgi=compose(f, g.segs[i]); + fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1])); + result.concat(fgi); + } + return result; +} + +template +Piecewise Piecewise::operator()(SBasis f){return compose((*this),f);} +template +Piecewise Piecewise::operator()(Piecewisef){return compose((*this),f);} + +template +Piecewise integral(Piecewise const &a) { + Piecewise 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 +Piecewise derivative(Piecewise const &a) { + Piecewise 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 roots(Piecewise 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 : diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp index 136687ae8..fdc1cefe5 100644 --- a/src/2geom/poly-dk-solve.cpp +++ b/src/2geom/poly-dk-solve.cpp @@ -1,64 +1,64 @@ -#include "poly-dk-solve.h" -#include - -/*** implementation of the Durand-Kerner method. seems buggy*/ - -std::complex evalu(Poly const & p, std::complex x) { - std::complex result = 0; - std::complex xx = 1; - - for(unsigned i = 0; i < p.size(); i++) { - result += p[i]*xx; - xx *= x; - } - return result; -} - -std::vector > DK(Poly const & ply, const double tol) { - std::vector > roots; - const int N = ply.degree(); - - std::complex b(0.4, 0.9); - std::complex 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 denom = 1; - std::complex 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 dr = evalu(ply, R)/denom; - error += norm(dr); - roots[r_i] = R - dr; - } - /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(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 : +#include "poly-dk-solve.h" +#include + +/*** implementation of the Durand-Kerner method. seems buggy*/ + +std::complex evalu(Poly const & p, std::complex x) { + std::complex result = 0; + std::complex xx = 1; + + for(unsigned i = 0; i < p.size(); i++) { + result += p[i]*xx; + xx *= x; + } + return result; +} + +std::vector > DK(Poly const & ply, const double tol) { + std::vector > roots; + const int N = ply.degree(); + + std::complex b(0.4, 0.9); + std::complex 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 denom = 1; + std::complex 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 dr = evalu(ply, R)/denom; + error += norm(dr); + roots[r_i] = R - dr; + } + /*std::copy(roots.begin(), roots.end(), std::ostream_iterator >(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 : diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp index fb3c2075b..766f16eda 100644 --- a/src/2geom/poly-laguerre-solve.cpp +++ b/src/2geom/poly-laguerre-solve.cpp @@ -1,147 +1,147 @@ -#include "poly-laguerre-solve.h" -#include - -typedef std::complex 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 -laguerre(Poly p, const double tol) { - std::vector 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 -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 : +#include "poly-laguerre-solve.h" +#include + +typedef std::complex 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 +laguerre(Poly p, const double tol) { + std::vector 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 +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 ef583c0b5..62c769c8b 100644 --- a/src/2geom/quadtree.h +++ b/src/2geom/quadtree.h @@ -1,4 +1,4 @@ -#include +#include #include class Quad{ @@ -20,7 +20,7 @@ public: 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 5bf0d2876..7157bc885 100644 --- a/src/2geom/sbasis.cpp +++ b/src/2geom/sbasis.cpp @@ -1,490 +1,490 @@ -/* - * sbasis.cpp - S-power basis function class + supporting classes - * - * Authors: - * Nathan Hurst - * Michael Sloan - * - * 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 - -#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; -} +/* + * sbasis.cpp - S-power basis function class + supporting classes + * + * Authors: + * Nathan Hurst + * Michael Sloan + * + * 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 + +#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; +} */ - -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= 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 << ": ---------" < 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= 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 << ": ---------" < - -namespace Geom { - -std::vector > sweep_bounds(std::vector rs) { - std::vector events; events.reserve(rs.size()*2); - std::vector > 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 open; - for(unsigned i = 0; i < events.size(); i++) { - unsigned ix = events[i].ix; - if(events[i].closing) { - std::vector::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 > sweep_bounds(std::vector a, std::vector b) { - std::vector > pairs(a.size()); - if(a.empty() || b.empty()) return pairs; - std::vector 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 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 > fake_cull(unsigned a, unsigned b) { - std::vector > ret; - - std::vector 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; -} - -} +#include "sweep.h" + +#include + +namespace Geom { + +std::vector > sweep_bounds(std::vector rs) { + std::vector events; events.reserve(rs.size()*2); + std::vector > 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 open; + for(unsigned i = 0; i < events.size(); i++) { + unsigned ix = events[i].ix; + if(events[i].closing) { + std::vector::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 > sweep_bounds(std::vector a, std::vector b) { + std::vector > pairs(a.size()); + if(a.empty() || b.empty()) return pairs; + std::vector 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 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 > fake_cull(unsigned a, unsigned b) { + std::vector > ret; + + std::vector 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; +} + +} diff --git a/src/live_effects/lpe-pathalongpath.cpp b/src/live_effects/lpe-pathalongpath.cpp index 4837317bf..9ac75cb3c 100644 --- a/src/live_effects/lpe-pathalongpath.cpp +++ b/src/live_effects/lpe-pathalongpath.cpp @@ -1,221 +1,221 @@ -#define INKSCAPE_LPE_PATHALONGPATH_CPP - -/* - * Copyright (C) Johan Engelen 2007 - * - * 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 -#include -#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 -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 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 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(&bend_path) ); -/* Delayed until 0.47 - registerParameter( dynamic_cast(&width_path) ); - registerParameter( dynamic_cast(&width_path_range) ); -*/ - registerParameter( dynamic_cast(©type) ); - registerParameter( dynamic_cast(&prop_scale) ); - registerParameter( dynamic_cast(&scale_y_rel) ); - registerParameter( dynamic_cast(&vertical_pattern) ); - - prop_scale.param_set_digits(3); - prop_scale.param_set_increments(0.01, 0.10); -} - -LPEPathAlongPath::~LPEPathAlongPath() -{ - -} - - -Geom::Piecewise > -LPEPathAlongPath::doEffect_pwd2 (Geom::Piecewise > & 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 > uskeleton = arc_length_parametrization(Piecewise >(bend_path),2,.1); - uskeleton = remove_short_cuts(uskeleton,.01); - Rect uskeletonbounds = bounds_exact(uskeleton); - uskeleton -= uskeletonbounds.midpoint(); - - Piecewise > n = rot90(derivative(uskeleton)); - n = force_continuity(remove_short_cuts(n,.1)); - - D2 > patternd2 = make_cuts_independant(pwd2_in); - Piecewise x = vertical_pattern.get_value() ? Piecewise(patternd2[1]) : Piecewise(patternd2[0]); - Piecewise y = vertical_pattern.get_value() ? Piecewise(patternd2[0]) : Piecewise(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 > widthpwd2 = arc_length_parametrization(Piecewise >(width_path),2,.1); - D2 > widthd2pw = make_cuts_independant(widthpwd2); - Piecewise width = (Piecewise(widthd2pw[Y]) - uskeletonbounds[Y].middle()) / width_path_range; -*/ - - double offs = 0; - Piecewise > output; - for (int i=0; i > pwd2; - std::vector 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 > 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( 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( 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 : +#define INKSCAPE_LPE_PATHALONGPATH_CPP + +/* + * Copyright (C) Johan Engelen 2007 + * + * 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 +#include +#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 +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 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 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(&bend_path) ); +/* Delayed until 0.47 + registerParameter( dynamic_cast(&width_path) ); + registerParameter( dynamic_cast(&width_path_range) ); +*/ + registerParameter( dynamic_cast(©type) ); + registerParameter( dynamic_cast(&prop_scale) ); + registerParameter( dynamic_cast(&scale_y_rel) ); + registerParameter( dynamic_cast(&vertical_pattern) ); + + prop_scale.param_set_digits(3); + prop_scale.param_set_increments(0.01, 0.10); +} + +LPEPathAlongPath::~LPEPathAlongPath() +{ + +} + + +Geom::Piecewise > +LPEPathAlongPath::doEffect_pwd2 (Geom::Piecewise > & 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 > uskeleton = arc_length_parametrization(Piecewise >(bend_path),2,.1); + uskeleton = remove_short_cuts(uskeleton,.01); + Rect uskeletonbounds = bounds_exact(uskeleton); + uskeleton -= uskeletonbounds.midpoint(); + + Piecewise > n = rot90(derivative(uskeleton)); + n = force_continuity(remove_short_cuts(n,.1)); + + D2 > patternd2 = make_cuts_independant(pwd2_in); + Piecewise x = vertical_pattern.get_value() ? Piecewise(patternd2[1]) : Piecewise(patternd2[0]); + Piecewise y = vertical_pattern.get_value() ? Piecewise(patternd2[0]) : Piecewise(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 > widthpwd2 = arc_length_parametrization(Piecewise >(width_path),2,.1); + D2 > widthd2pw = make_cuts_independant(widthpwd2); + Piecewise width = (Piecewise(widthd2pw[Y]) - uskeletonbounds[Y].middle()) / width_path_range; +*/ + + double offs = 0; + Piecewise > output; + for (int i=0; i > pwd2; + std::vector 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 > 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( 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( 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 : diff --git a/src/live_effects/lpe-pathalongpath.h b/src/live_effects/lpe-pathalongpath.h index be064dd04..b556efd88 100644 --- a/src/live_effects/lpe-pathalongpath.h +++ b/src/live_effects/lpe-pathalongpath.h @@ -1,55 +1,55 @@ -#ifndef INKSCAPE_LPE_PATHALONGPATH_H -#define INKSCAPE_LPE_PATHALONGPATH_H - -/* - * Inkscape::LPEPathAlongPath - * -* Copyright (C) Johan Engelen 2007 - * - * 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 > doEffect_pwd2 (Geom::Piecewise > & pwd2_in); - - virtual void resetDefaults(SPItem * item); - -private: - PathParam bend_path; - //PathParam width_path; - //ScalarParam width_path_range; - EnumParam 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 +#ifndef INKSCAPE_LPE_PATHALONGPATH_H +#define INKSCAPE_LPE_PATHALONGPATH_H + +/* + * Inkscape::LPEPathAlongPath + * +* Copyright (C) Johan Engelen 2007 + * + * 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 > doEffect_pwd2 (Geom::Piecewise > & pwd2_in); + + virtual void resetDefaults(SPItem * item); + +private: + PathParam bend_path; + //PathParam width_path; + //ScalarParam width_path_range; + EnumParam 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 -- 2.30.2