src/2geom/piecewise.cpp

   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  */
  31
  32 #include <2geom/piecewise.h>
  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]);
 161
 162         for (unsigned r=0; r<rts.size(); r++){
 163             result.push_back(f.mapToDomain(rts[r], i));
 164         }
 165     }
 166     return result;
 167 }
 168
 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 }
 181
 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 :