src/2geom/d2-sbasis.cpp

   1 #include <2geom/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_independent(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 /** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
  62  * @return  \f[
  63  *      f(t) \rightarrow \left\{
  64  *      \begin{array}{c}
  65  *      a_1 \bullet b_1 \\
  66  *      a_2 \bullet b_2 \\
  67  *      \ldots \\
  68  *      a_n \bullet b_n \\
  69  *      \end{array}\right.
  70  * \f]
  71  * @relates Piecewise */
  72 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b)
  73 {
  74   Piecewise<SBasis > result;
  75   if (a.empty() || b.empty()) return result;
  76   Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
  77   Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
  78
  79   result.push_cut(aa.cuts.front());
  80   for (unsigned i=0; i<aa.size(); i++){
  81     result.push(dot(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
  82   }
  83   return result;
  84 }
  85
  86 /** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
  87  * @return  \f[
  88  *      f(t) \rightarrow \left\{
  89  *      \begin{array}{c}
  90  *      a_1 \bullet b \\
  91  *      a_2 \bullet b \\
  92  *      \ldots \\
  93  *      a_n \bullet b \\
  94  *      \end{array}\right.
  95  * \f]
  96  * @relates Piecewise */
  97 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Point const &b)
  98 {
  99   Piecewise<SBasis > result;
 100   if (a.empty()) return result;
 101
 102   result.push_cut(a.cuts.front());
 103   for (unsigned i = 0; i < a.size(); ++i){
 104     result.push(dot(a.segs[i],b), a.cuts[i+1]);
 105   }
 106   return result;
 107 }
 108
 109
 110 Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a,
 111                         Piecewise<D2<SBasis> > const &b){
 112   Piecewise<SBasis > result;
 113   if (a.empty() || b.empty()) return result;
 114   Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
 115   Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
 116
 117   result.push_cut(aa.cuts.front());
 118   for (unsigned i=0; i<a.size(); i++){
 119     result.push(cross(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
 120   }
 121   return result;
 122 }
 123
 124 Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Matrix const &m) {
 125   Piecewise<D2<SBasis> > result;
 126   if(a.empty()) return result;
 127   result.push_cut(a.cuts[0]);
 128   for (unsigned i = 0; i < a.size(); i++) {
 129     result.push(a[i] * m, a.cuts[i+1]);
 130   }
 131   return result;
 132 }
 133
 134 //if tol>0, only force continuity where the jump is smaller than tol.
 135 Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double tol, bool closed)
 136 {
 137     if (f.size()==0) return f;
 138     Piecewise<D2<SBasis> > result=f;
 139     unsigned cur   = (closed)? 0:1;
 140     unsigned prev  = (closed)? f.size()-1:0;
 141     while(cur<f.size()){
 142         Point pt0 = f.segs[prev].at1();
 143         Point pt1 = f.segs[cur ].at0();
 144         if (tol<=0 || L2sq(pt0-pt1)<tol*tol){
 145             pt0 = (pt0+pt1)/2;
 146             for (unsigned dim=0; dim<2; dim++){
 147                 SBasis &prev_sb=result.segs[prev][dim];
 148                 SBasis &cur_sb =result.segs[cur][dim];
 149                 Coord const c=pt0[dim];
 150                 if (prev_sb.empty()) {
 151                   prev_sb = SBasis(Linear(0.0, c));
 152                 } else {
 153                   prev_sb[0][1] = c;
 154                 }
 155                 if (cur_sb.empty()) {
 156                   cur_sb = SBasis(Linear(c, 0.0));
 157                 } else {
 158                   cur_sb[0][0] = c;
 159                 }
 160             }
 161         }
 162         prev = cur++;
 163     }
 164     return result;
 165 }
 166
 167 std::vector<Geom::Piecewise<Geom::D2<Geom::SBasis> > >
 168 split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbin, double tol)
 169 {
 170     using namespace Geom;
 171     std::vector<Piecewise<D2<SBasis> > > ret;
 172     unsigned piece_start = 0;
 173     for (unsigned i=0; i<pwsbin.segs.size(); i++){
 174         if (i==(pwsbin.segs.size()-1) || L2(pwsbin.segs[i].at1()- pwsbin.segs[i+1].at0()) > tol){
 175             Piecewise<D2<SBasis> > piece;
 176             piece.cuts.push_back(pwsbin.cuts[piece_start]);
 177             for (unsigned j = piece_start; j<i+1; j++){
 178                 piece.segs.push_back(pwsbin.segs[j]);
 179                 piece.cuts.push_back(pwsbin.cuts[j+1]);
 180             }
 181             ret.push_back(piece);
 182             piece_start = i+1;
 183         }
 184     }
 185     return ret;
 186 }
 187
 188 static void set_first_point(Piecewise<D2<SBasis> > &f, Point a){
 189     if ( f.empty() ){
 190         f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(Linear(a[X]),Linear(a[Y]))));
 191         return;
 192     }
 193     for (unsigned dim=0; dim<2; dim++){
 194         if (f.segs.front()[dim].size() == 0){
 195             f.segs.front()[dim] = SBasis(Linear(a[dim],0));
 196         }else{
 197             f.segs.front()[dim][0][0] = a[dim];
 198         }
 199     }
 200 }
 201 static void set_last_point(Piecewise<D2<SBasis> > &f, Point a){
 202     if ( f.empty() ){
 203         f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(Linear(a[X]),Linear(a[Y]))));
 204         return;
 205     }
 206     for (unsigned dim=0; dim<2; dim++){
 207         if (f.segs.back()[dim].size() == 0){
 208             f.segs.back()[dim] = SBasis(Linear(0,a[dim]));
 209         }else{
 210             f.segs.back()[dim][0][1] = a[dim];
 211         }
 212     }
 213 }
 214
 215 std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol){
 216
 217     if ( f.size()==0 ) return f;
 218     std::vector<Piecewise<D2<SBasis> > > result;
 219     std::vector<std::vector<unsigned> > pre_result;
 220     for (unsigned i=0; i<f.size(); i++){
 221         bool inserted = false;
 222         Point a = f[i].firstValue();
 223         Point b = f[i].lastValue();
 224         for (unsigned j=0; j<pre_result.size(); j++){
 225             Point aj = f.at(pre_result[j].back()).lastValue();
 226             Point bj = f.at(pre_result[j].front()).firstValue();
 227             if ( L2(a-aj) < tol ) {
 228                 pre_result[j].push_back(i);
 229                 inserted = true;
 230                 break;
 231             }
 232             if ( L2(b-bj) < tol ) {
 233                 pre_result[j].insert(pre_result[j].begin(),i);
 234                 inserted = true;
 235                 break;
 236             }
 237         }
 238         if (!inserted) {
 239             pre_result.push_back(std::vector<unsigned>());
 240             pre_result.back().push_back(i);
 241         }
 242     }
 243     for (unsigned i=0; i<pre_result.size(); i++){
 244         Piecewise<D2<SBasis> > comp;
 245         for (unsigned j=0; j<pre_result[i].size(); j++){
 246             Piecewise<D2<SBasis> > new_comp = f.at(pre_result[i][j]);
 247             if ( j>0 ){
 248                 set_first_point( new_comp, comp.segs.back().at1() );
 249             }
 250             comp.concat(new_comp);
 251         }
 252         if ( L2(comp.firstValue()-comp.lastValue()) < tol ){
 253             //TODO: check sizes!!!
 254             set_last_point( comp, comp.segs.front().at0() );
 255         }
 256         result.push_back(comp);
 257     }
 258     return result;
 259     return f;
 260 }
 261
 262
 263 }  // namespace Geom
 264
 265
 266 /*
 267   Local Variables:
 268   mode:c++
 269   c-file-style:"stroustrup"
 270   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 271   indent-tabs-mode:nil
 272   fill-column:99
 273   End:
 274 */
 275 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :