src/2geom/d2-sbasis.cpp

   1 #include "d2.h"
   2 /* One would think that we would include d2-sbasis.h, however,
   3  * you cannot actually include it in anything - only d2 may import it.
   4  * This is due to the trickinesses of template submatching. */
   5
   6 namespace Geom {
   7
   8 SBasis L2(D2<SBasis> const & a, unsigned k) { return sqrt(dot(a, a), k); }
   9
  10 D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b) {
  11     return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));
  12 }
  13
  14 D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b) {
  15     return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));
  16 }
  17
  18 D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms) {
  19     return D2<SBasis>(truncate(a[X], terms), truncate(a[Y], terms));
  20 }
  21
  22 unsigned sbasis_size(D2<SBasis> const & a) {
  23     return std::max((unsigned) a[0].size(), (unsigned) a[1].size());
  24 }
  25
  26 //TODO: Is this sensical? shouldn't it be like pythagorean or something?
  27 double tail_error(D2<SBasis> const & a, unsigned tail) {
  28     return std::max(a[0].tailError(tail), a[1].tailError(tail));
  29 }
  30
  31 Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a) {
  32     Piecewise<SBasis> x = partition(a[0], a[1].cuts), y = partition(a[1], a[0].cuts);
  33     assert(x.size() == y.size());
  34     Piecewise<D2<SBasis> > ret;
  35     for(unsigned i = 0; i < x.size(); i++)
  36         ret.push_seg(D2<SBasis>(x[i], y[i]));
  37     ret.cuts.insert(ret.cuts.end(), x.cuts.begin(), x.cuts.end());
  38     return ret;
  39 }
  40
  41 D2<Piecewise<SBasis> > make_cuts_independant(Piecewise<D2<SBasis> > const &a) {
  42     D2<Piecewise<SBasis> > ret;
  43     for(unsigned d = 0; d < 2; d++) {
  44         for(unsigned i = 0; i < a.size(); i++)
  45             ret[d].push_seg(a[i][d]);
  46         ret[d].cuts.insert(ret[d].cuts.end(), a.cuts.begin(), a.cuts.end());
  47     }
  48     return ret;
  49 }
  50
  51 Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &M){
  52   Piecewise<D2<SBasis> > result;
  53   if (M.empty()) return M;
  54   result.push_cut(M.cuts[0]);
  55   for (unsigned i=0; i<M.size(); i++){
  56     result.push(rot90(M[i]),M.cuts[i+1]);
  57   }
  58   return result;
  59 }
  60
  61 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a,
  62                       Piecewise<D2<SBasis> > const &b){
  63   Piecewise<SBasis > result;
  64   if (a.empty() || b.empty()) return result;
  65   Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
  66   Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
  67
  68   result.push_cut(aa.cuts.front());
  69   for (unsigned i=0; i<aa.size(); i++){
  70     result.push(dot(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
  71   }
  72   return result;
  73 }
  74
  75 Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a,
  76                         Piecewise<D2<SBasis> > const &b){
  77   Piecewise<SBasis > result;
  78   if (a.empty() || b.empty()) return result;
  79   Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
  80   Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
  81
  82   result.push_cut(aa.cuts.front());
  83   for (unsigned i=0; i<a.size(); i++){
  84     result.push(cross(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
  85   }
  86   return result;
  87 }
  88
  89 Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Matrix const &m) {
  90   Piecewise<D2<SBasis> > result;
  91   if(a.empty()) return result;
  92   result.push_cut(a.cuts[0]);
  93   for (unsigned i = 0; i < a.size(); i++) {
  94     result.push(a[i] * m, a.cuts[i+1]);
  95   }
  96   return result;
  97 }
  98
  99 /* Replaced by remove_short_cuts in piecewise.h
 100 //this recursively removes the shortest cut interval until none is shorter than tol.
 101 //TODO: code this in a more efficient way!
 102 Piecewise<D2<SBasis> > remove_short_cuts(Piecewise<D2<SBasis> > const &f, double tol){
 103     double min = tol;
 104     unsigned idx = f.size();
 105     for(unsigned i=0; i<f.size(); i++){
 106         if (min > f.cuts[i+1]-f.cuts[i]){
 107             min = f.cuts[i+1]-f.cuts[i];
 108             idx = int(i);
 109         }
 110     }
 111     if (idx==f.size()){
 112         return f;
 113     }
 114     if (f.size()==1) {
 115         //removing this seg would result in an empty pw<d2<sb>>...
 116         return f;
 117     }
 118     Piecewise<D2<SBasis> > new_f=f;
 119     for (int dim=0; dim<2; dim++){
 120         double v = Hat(f.segs.at(idx)[dim][0]);
 121         //TODO: what about closed curves?
 122         if (idx>0 && f.segs.at(idx-1).at1()==f.segs.at(idx).at0())
 123             new_f.segs.at(idx-1)[dim][0][1] = v;
 124         if (idx<f.size() && f.segs.at(idx+1).at0()==f.segs.at(idx).at1())
 125             new_f.segs.at(idx+1)[dim][0][0] = v;
 126     }
 127     double t = (f.cuts.at(idx)+f.cuts.at(idx+1))/2;
 128     new_f.cuts.at(idx+1) = t;
 129
 130     new_f.segs.erase(new_f.segs.begin()+idx);
 131     new_f.cuts.erase(new_f.cuts.begin()+idx);
 132     return remove_short_cuts(new_f, tol);
 133 }
 134 */
 135
 136 //if tol>0, only force continuity where the jump is smaller than tol.
 137 Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f,
 138                                         double tol,
 139                                         bool closed){
 140     if (f.size()==0) return f;
 141     Piecewise<D2<SBasis> > result=f;
 142     unsigned cur   = (closed)? 0:1;
 143     unsigned prev  = (closed)? f.size()-1:0;
 144     while(cur<f.size()){
 145         Point pt0 = f.segs[prev].at1();
 146         Point pt1 = f.segs[cur ].at0();
 147         if (tol<=0 || L2sq(pt0-pt1)<tol*tol){
 148             pt0 = (pt0+pt1)/2;
 149             for (unsigned dim=0; dim<2; dim++){
 150                 SBasis &prev_sb=result.segs[prev][dim];
 151                 SBasis &cur_sb =result.segs[cur][dim];
 152                 Coord const c=pt0[dim];
 153                 if (prev_sb.empty()) {
 154                   prev_sb.push_back(Linear(0.0, c));
 155                 } else {
 156                   prev_sb[0][1] = c;
 157                 }
 158                 if (cur_sb.empty()) {
 159                   cur_sb.push_back(Linear(c, 0.0));
 160                 } else {
 161                   cur_sb[0][0] = c;
 162                 }
 163             }
 164         }
 165         prev = cur++;
 166     }
 167     return result;
 168 }
 169 }