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
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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 */
32 #include <2geom/piecewise.h>
33 #include <iterator>
34 #include <map>
36 namespace Geom {
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 }
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 }
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)%...
77 k+=1;
78 r.resize(k, Linear(0,0));
79 c.resize(k, Linear(0,0));
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 }
88 if (r.tailError(k)<tol) return Piecewise<SBasis>(c);
89 }
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 }
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 */
107 std::map<double,unsigned> compose_pullback(std::vector<double> const &values, SBasis const &g){
108 std::map<double,unsigned> result;
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 }
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 }
152 //move idx back from levels f.cuts
153 idx+=1;
154 return idx;
155 }
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]);
162 for (unsigned r=0; r<rts.size(); r++){
163 result.push_back(f.mapToDomain(rts[r], i));
164 }
165 }
166 return result;
167 }
169 std::vector<std::vector<double> > multi_roots(Piecewise<SBasis> const &f, std::vector<double> const &values) {
170 std::vector<std::vector<double> > result(values.size());
171 for (unsigned i=0; i<f.size(); i++) {
172 std::vector<std::vector<double> > rts = multi_roots(f.segs[i], values);
173 for(unsigned j=0; j<rts.size(); j++) {
174 for(unsigned r=0; r<rts[j].size(); r++){
175 result[j].push_back(f.mapToDomain(rts[j][r], i));
176 }
177 }
178 }
179 return result;
180 }
182 }
183 /*
184 Local Variables:
185 mode:c++
186 c-file-style:"stroustrup"
187 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
188 indent-tabs-mode:nil
189 fill-column:99
190 End:
191 */
192 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :