Rev 15347 | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
12065 | fschmid | 1 | /* |
2 | * piecewise.cpp - Piecewise function class |
||
3 | * |
||
4 | * Copyright 2007 Michael Sloan <mgsloan@gmail.com> |
||
5 | * Copyright 2007 JF Barraud |
||
6 | * |
||
7 | * This library is free software; you can redistribute it and/or |
||
8 | * modify it either under the terms of the GNU Lesser General Public |
||
9 | * License version 2.1 as published by the Free Software Foundation |
||
10 | * (the "LGPL") or, at your option, under the terms of the Mozilla |
||
11 | * Public License Version 1.1 (the "MPL"). If you do not alter this |
||
12 | * notice, a recipient may use your version of this file under either |
||
13 | * the MPL or the LGPL. |
||
14 | * |
||
15 | * You should have received a copy of the LGPL along with this library |
||
16 | * in the file COPYING-LGPL-2.1; if not, output to the Free Software |
||
18122 | mrdocs | 17 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
12065 | fschmid | 18 | * You should have received a copy of the MPL along with this library |
19 | * in the file COPYING-MPL-1.1 |
||
20 | * |
||
21 | * The contents of this file are subject to the Mozilla Public License |
||
22 | * Version 1.1 (the "License"); you may not use this file except in |
||
23 | * compliance with the License. You may obtain a copy of the License at |
||
24 | * http://www.mozilla.org/MPL/ |
||
25 | * |
||
26 | * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY |
||
27 | * OF ANY KIND, either express or implied. See the LGPL or the MPL for |
||
28 | * the specific language governing rights and limitations. |
||
29 | * |
||
30 | */ |
||
31 | |||
15347 | fschmid | 32 | #include "piecewise.h" |
12065 | fschmid | 33 | #include <iterator> |
34 | #include <map> |
||
35 | |||
36 | namespace Geom { |
||
37 | |||
38 | Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k) { |
||
39 | Piecewise<SBasis> pa = partition(a, b.cuts), pb = partition(b, a.cuts); |
||
40 | Piecewise<SBasis> ret = Piecewise<SBasis>(); |
||
41 | assert(pa.size() == pb.size()); |
||
42 | ret.cuts = pa.cuts; |
||
43 | for (unsigned i = 0; i < pa.size(); i++) |
||
44 | ret.push_seg(divide(pa[i], pb[i], k)); |
||
45 | return ret; |
||
46 | } |
||
47 | |||
48 | Piecewise<SBasis> |
||
49 | divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero) { |
||
50 | Piecewise<SBasis> pa = partition(a, b.cuts), pb = partition(b, a.cuts); |
||
51 | Piecewise<SBasis> ret = Piecewise<SBasis>(); |
||
52 | assert(pa.size() == pb.size()); |
||
53 | for (unsigned i = 0; i < pa.size(); i++){ |
||
54 | Piecewise<SBasis> divi = divide(pa[i], pb[i], tol, k, zero); |
||
55 | divi.setDomain(Interval(pa.cuts[i],pa.cuts[i+1])); |
||
56 | ret.concat(divi); |
||
57 | } |
||
58 | return ret; |
||
59 | } |
||
60 | Piecewise<SBasis> divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero){ |
||
61 | return divide(a,Piecewise<SBasis>(b),tol,k,zero); |
||
62 | } |
||
63 | Piecewise<SBasis> divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero){ |
||
64 | return divide(Piecewise<SBasis>(a),b,tol,k,zero); |
||
65 | } |
||
66 | Piecewise<SBasis> divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero) { |
||
67 | if (b.tailError(0)<2*zero){ |
||
68 | //TODO: have a better look at sgn(b). |
||
69 | double sgn= (b(.5)<0.)?-1.:1; |
||
70 | return Piecewise<SBasis>(Linear(sgn/zero)*a); |
||
71 | } |
||
72 | |||
73 | if (fabs(b.at0())>zero && fabs(b.at1())>zero ){ |
||
74 | SBasis c,r=a; |
||
75 | //TODO: what is a good relative tol? atm, c=a/b +/- (tol/a)%... |
||
76 | |||
77 | k+=1; |
||
78 | r.resize(k, Linear(0,0)); |
||
79 | c.resize(k, Linear(0,0)); |
||
80 | |||
81 | //assert(b.at0()!=0 && b.at1()!=0); |
||
82 | for (unsigned i=0; i<k; i++){ |
||
83 | Linear ci = Linear(r[i][0]/b[0][0],r[i][1]/b[0][1]); |
||
84 | c[i]=ci; |
||
85 | r-=shift(ci*b,i); |
||
86 | } |
||
87 | |||
88 | if (r.tailError(k)<tol) return Piecewise<SBasis>(c); |
||
89 | } |
||
90 | |||
91 | Piecewise<SBasis> c0,c1; |
||
92 | c0 = divide(compose(a,Linear(0.,.5)),compose(b,Linear(0.,.5)),tol,k); |
||
93 | c1 = divide(compose(a,Linear(.5,1.)),compose(b,Linear(.5,1.)),tol,k); |
||
94 | c0.setDomain(Interval(0.,.5)); |
||
95 | c1.setDomain(Interval(.5,1.)); |
||
96 | c0.concat(c1); |
||
97 | return c0; |
||
98 | } |
||
99 | |||
100 | |||
101 | //-- compose(pw<T>,SBasis) --------------- |
||
102 | /* |
||
103 | the purpose of the following functions is only to reduce the code in piecewise.h |
||
104 | TODO: use a vector<pairs<double,unsigned> > instead of a map<double,unsigned>. |
||
105 | */ |
||
106 | |||
107 | std::map<double,unsigned> compose_pullback(std::vector<double> const &values, SBasis const &g){ |
||
108 | std::map<double,unsigned> result; |
||
109 | |||
110 | std::vector<std::vector<double> > roots = multi_roots(g, values); |
||
111 | for(unsigned i=0; i<roots.size(); i++){ |
||
112 | for(unsigned j=0; j<roots[i].size();j++){ |
||
113 | result[roots[i][j]]=i; |
||
114 | } |
||
115 | } |
||
116 | // Also map 0 and 1 to the first value above(or =) g(0) and g(1). |
||
117 | if(result.count(0.)==0){ |
||
118 | unsigned i=0; |
||
119 | while (i<values.size()&&(g.at0()>values[i])) i++; |
||
120 | result[0.]=i; |
||
121 | } |
||
122 | if(result.count(1.)==0){ |
||
123 | unsigned i=0; |
||
124 | while (i<values.size()&&(g.at1()>values[i])) i++; |
||
125 | result[1.]=i; |
||
126 | } |
||
127 | return(result); |
||
128 | } |
||
129 | |||
130 | int compose_findSegIdx(std::map<double,unsigned>::iterator const &cut, |
||
131 | std::map<double,unsigned>::iterator const &next, |
||
132 | std::vector<double> const &levels, |
||
133 | SBasis const &g){ |
||
134 | double t0=(*cut).first; |
||
135 | unsigned idx0=(*cut).second; |
||
136 | double t1=(*next).first; |
||
137 | unsigned idx1=(*next).second; |
||
138 | assert(t0<t1); |
||
139 | int idx; //idx of the relevant f.segs |
||
140 | if (std::max(idx0,idx1)==levels.size()){ //g([t0,t1]) is above the top level, |
||
141 | idx=levels.size()-1; |
||
142 | } else if (idx0 != idx1){ //g([t0,t1]) crosses from level idx0 to idx1, |
||
143 | idx=std::min(idx0,idx1); |
||
144 | } else if(g((t0+t1)/2) < levels[idx0]) { //g([t0,t1]) is a 'U' under level idx0, |
||
145 | idx=idx0-1; |
||
146 | } else if(g((t0+t1)/2) > levels[idx0]) { //g([t0,t1]) is a 'bump' over level idx0, |
||
147 | idx=idx0; |
||
148 | } else { //g([t0,t1]) is contained in level idx0!... |
||
149 | idx = (idx0==levels.size())? idx0-1:idx0; |
||
150 | } |
||
151 | |||
152 | //move idx back from levels f.cuts |
||
153 | idx+=1; |
||
154 | return idx; |
||
155 | } |
||
156 | |||
157 | std::vector<double> roots(Piecewise<SBasis> const &f){ |
||
158 | std::vector<double> result; |
||
159 | for (unsigned i=0; i<f.size(); i++){ |
||
160 | std::vector<double> rts=roots(f.segs[i]); |
||
15347 | fschmid | 161 | rts=roots(f.segs[i]); |
12065 | fschmid | 162 | |
163 | for (unsigned r=0; r<rts.size(); r++){ |
||
164 | result.push_back(f.mapToDomain(rts[r], i)); |
||
165 | } |
||
166 | } |
||
167 | return result; |
||
168 | } |
||
169 | |||
170 | } |
||
171 | /* |
||
172 | Local Variables: |
||
173 | mode:c++ |
||
174 | c-file-style:"stroustrup" |
||
175 | c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) |
||
176 | indent-tabs-mode:nil |
||
177 | fill-column:99 |
||
178 | End: |
||
179 | */ |
||
180 | // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : |