---

TensorBlossom.H

Go to the documentation of this file.

#ifndef _TENSORBLOSSOM_H
#define _TENSORBLOSSOM_H

#include <xstl.H>
#include <TensorManifold.H>
#include <Vector.H>
#include <Combinatory.H>

namespace xchen
{
  
  enum { ruled, revolution, sweep };                              

  template <int dims, int sp_dims>
  struct generate_iso_functor;
  template <int dims, int sp_dims>
  struct draw_iso_functor;

  template<int dims=2, int sp_dims=3>
  class TensorBlossom : public TensorManifold<dims, sp_dims>
  {
      template<int dims2, int sp_dims2> friend class TensorBlossom;
      friend class glModel;
      template <int dims2, int sp_dims2> friend struct generate_iso_functor;
      template <int dims2, int sp_dims2> friend struct draw_iso_functor;

      def_type_name_for_tensor_manifold(dims, sp_dims);
      typedef TensorBlossom<dims, sp_dims> MyType;
      typedef TensorManifold<dims, sp_dims> BaseType;
      typedef TensorBlossom<dims-1,sp_dims> IsoMFType;

#define SubMFType TensorBlossom<dims-1,from_sp_dims>

  private:
      knot_t ctl_kv;                                              // defining data: ctl_kv is the defining knot vectors.
      bool periodic[dims];                                        // defining data: periodic end condtions.

      vector<double> dmn[dims];                                   // derived data: break knots extracted from ctl_kv.
      bool domain_computed;                                       // For now assume unmutated after computed.
      
      knot_t sbd_kv;                                              // subdivision data: sbd_kv is for knot insertion and subdvion algo.
      vector< IsoMFType > iso_mf[dims];                           // iso sub_manifold at break points in each direction.
      
  public:
      template<int from_sp_dims>
      TensorBlossom(int mtd,SubMFType const*,SubMFType const* =0);

      TensorBlossom(int arg1, ...);                               
      TensorBlossom(mesh_t const& mesh);                          
      TensorBlossom& Transform(xchen::Transform const& trans);    
      
      knot_t& KnotVector() { return ctl_kv; }                            
      bool Insert(double u,int i=0);                                     
      bool IsPeriodicEnd(int i=0)           { return periodic[i]; }      
      void SetPeriodicEnd(int i=0,bool v=1) {periodic[i]=v; setKt(); }   
      void SetUniformFloatingEndKnots();                                 
      void SetUniformOpenEndKnots();                                     

      IsoMFType GetSubTensorBlossom(double u, int=0);                    

      void RaiseDegree(int inc=1, int i=0);                              
      void DecVarity(double*u){ForEachIofDims DecVarity(u[i],i);}        
      void DecVarity(double u, int i=0);                                 

      void Convert2FloatingEnd(int i)  { periodic2floating(i, ctldata);} 
      void Convert2Bzr(int i)          { convert2Bzr(i, ctldata); }      
      void Convert2OpenEnd(int i)      { convert2OpenEnd(i, ctldata); }  
      void Convert2Bzr()               { convert2Bzr(ctldata); }         
      void Convert2OpenEnd()           { convert2OpenEnd(ctldata); }     
      void Convert2FloatingEnd()       { periodic2floating(ctldata); }   

      void Subdivide(int depth=1)const;                                  

      friend ostream& operator<< <>(ostream& os, TensorBlossom const& rhs);

  private:
      template<int from_sp_dims> void construct_revolution_manifold(SubMFType const& rev_from);
      template<int from_sp_dims> void construct_ruled_manifold(SubMFType const& mf1, SubMFType const& mf2);

      enum { sbddata, ctldata };
      mesh_t& cur_msh(int applyon)                          { return applyon==sbddata? sbd_msh : ctl_msh; }
      knot_t& cur_kv(int applyon)                           { return applyon==sbddata? sbd_kv : ctl_kv; }

      bool insert(double u, int i, int applyon);      
      void blend(IdxRange const&, double u, int applyon);    // affine interpolation with u value at the given range.
      bool is_valid_range(IdxRange const& range, int on);
      
      void setKt() { ctl_kv.Clear(); ForEachIofDims ctl_kv[i].resize(periodic[i]? ctl_msh.GetSize(i)+1 : ctl_msh.GetSize(i) + deg[i] - 1); set_dflt_kt(); }
      void set_default_degs()                               { ForEachIofDims deg[i] = ctl_msh.GetSize(i)-1; }

      void init_sbd();                                      // init sbd data, and convert them to composite bzr.
      void init_sbd_knotvector(int i);                      // copy from ctl_kv and change to floating if periodic.
      void compact_sbd_kv();                                // sbd_kv is bzr form, delete the extra values.
      void subdivide_kv();                                  // insert middle value.
      void computeSegmentDomains();                         // Compute all break points for each direction based on ctl_kotvector, degree.
      void clean_sbd()                                      { BaseType::clean_sbd(); sbd_kv.Clear(); }

      void compute_1st_order_derivative();                  // all 1st order partial derivatives.
      void identifyEndPntsForPeriodicEnd(int i,int applyon);// periodic is implemented as open end. so this is to make sure no artificial gap.

      void periodic2floating(int i, int applyon);
      void periodic2floating(int applyon)                   { ForEachIofDims periodic2floating(i, applyon);  }       
      void convert2Bzr(int i, int applyon);                 // changes only made to subdivision data.
      void convert2Bzr(int applyon)                         { ForEachIofDims convert2Bzr(i, applyon);  }       
      void convert2OpenEnd(int i, int applyon);             // changes only made to subdivision data.
      void convert2OpenEnd(int applyon)                     { ForEachIofDims convert2OpenEnd(i, applyon);  }       

      void popFront(int i, int on)                          { cur_msh(on).PopFront(i); cur_kv(on)[i].erase(cur_kv(on)[i].begin());}//Erase front pnt&knot at i.
      void popBack(int i, int on)                           { cur_msh(on).PopBack(i); cur_kv(on)[i].pop_back(); }                  //Erase back pnt&knot at i.

      void unset_periodic()                                 { bzero(periodic, sizeof(bool) * dims); }
      void set_dflt_kt()                                    { SetUniformFloatingEndKnots(); }

      TensorBlossom(int deg[dims], int sz[dims]); 
      void GenerateUI()                               { BaseType::GenerateUI(); ui->AddboolControl("Draw Iso SubManifolds", 'i'); }   
      void Draw() const                               { BaseType::Draw();if(ui->GetboolControl("Draw Iso SubManifolds"))((MyType*)this)->draw_iso(); }

      void draw_iso();
      void generate_iso();
      void init_iso_generator()                             { iso_generator.blsm = this; }
      generate_iso_functor<dims, sp_dims> iso_generator;
  };


  template<int sp_dims> class TensorBlossom<0,sp_dims> { };

  
  template<int dims, int sp_dims>
  struct generate_iso_functor
  {
      bool first;    
      TensorBlossom<dims, sp_dims>* blsm;
      generate_iso_functor(TensorBlossom<dims, sp_dims>* _blsm=0) : first(true), blsm(_blsm) { }
      void operator()()
      {
        blsm->init_sbd();
        for(int i=0; i<dims; i++)
        {
          int sz[dims], dg[dims];
          blsm->sbd_msh.GetSize(sz);
          copy(sz+i+1, sz+dims, sz+i);
          copy(&blsm->deg[0], &blsm->deg[0]+i, dg);
          copy(&blsm->deg[0]+i+1, &blsm->deg[0]+dims, dg);

          // draw an iso every 2 * deg bzr segment. skip those already drawn by the previous calls.
          typename TensorBlossom<dims, sp_dims>::slc_iterator itr = blsm->sbd_msh.SliceIterate(i) + 2 * blsm->deg[i] * (first? 0 : 1); 
          int stride = (first? 2 : 4) * blsm->deg[i];
          for(; !itr.PastEnd(); itr += stride )
          {
            blsm->iso_mf[i].push_back( TensorBlossom<dims-1, sp_dims>(dg, sz) );
            copy(itr.Begin(), itr.End(), blsm->iso_mf[i].back().ControlMesh().Begin());    // copy the slice as the control mesh of the iso.
            blsm->iso_mf[i].back().rational = blsm->rational;
            
            int k = 0;
            for(int j=0; j<dims; j++)
            {
              if(j==i) continue;
              vector<double>::iterator out_itr = blsm->iso_mf[i].back().KnotVector()[k++].begin()-1;
              for(vector<double>::iterator itr = blsm->sbd_kv[j].begin(); itr != blsm->sbd_kv[j].end(); ++itr)
              {
                for(int i=0; i<blsm->deg[i]; i++)
                  *(++out_itr) = *itr;        // expand the bzr compact form of sbd_kv.
              }
            }
          }
          for(typename vector< TensorBlossom<dims-1, sp_dims> >::iterator itr = blsm->iso_mf[i].begin(); itr != blsm->iso_mf[i].end(); ++itr)
          {
            itr->Subdivide();
          }
        }
        first = false;
      }
  };
  
  template<int dims, int sp_dims>
  struct draw_iso_functor
  {
      TensorBlossom<dims, sp_dims>& blsm;
      draw_iso_functor(TensorBlossom<dims, sp_dims>& _blsm) : blsm(_blsm) { }
      void operator()()
      {
        glPushAttrib(GL_LIGHTING_BIT);
        glDisable(GL_LIGHTING);
        wglShading(RGBAS(1.0,1.0,1.0,1.0,1.0));
        for(int i=0; i<dims; i++)
        {
          for(typename vector< TensorBlossom<dims-1, sp_dims> > :: iterator itr = blsm.iso_mf[i].begin(); itr != blsm.iso_mf[i].end(); ++itr)
          {
            itr->draw(itr->sbd_msh);
          }
        }
        glPopAttrib();
      }
  };
  
  
  template<int sp_dims>
  struct generate_iso_functor<1, sp_dims>
  {
      TensorBlossom<1, sp_dims>* blsm;
      generate_iso_functor(TensorBlossom<1, sp_dims>* =0) { }
      void operator()() { }
  };
  template<int sp_dims>
  struct draw_iso_functor<1, sp_dims>
  {
      draw_iso_functor(TensorBlossom<1, sp_dims>&) { }
      void operator()() { }
  };
  
  

  
  template<int dims, int sp_dims>
  inline TensorBlossom<dims, sp_dims> :: TensorBlossom(mesh_t const& m) : BaseType(m), domain_computed(false)
  {
    set_default_degs();
    unset_periodic();
    setKt();
    SetUniformOpenEndKnots();
    set2bzr_srf_type();
    init_iso_generator();
  }
  template<int dims, int sp_dims>
  inline TensorBlossom<dims, sp_dims> :: TensorBlossom(int arg1, ...) : domain_computed(false) 
  { 
    int args[2*dims];
    int sz[dims];
    get_reverse_ints_va(args, 2*dims); 
    for(int i=dims-1; i>=0; i--)
    {
      deg[i] = args[2*i+1];
      sz[i] = args[2*i];
    }
    ctl_msh.Resize(sz);

    unset_periodic();
    setKt();
    set2bsp_srf_type();
    init_iso_generator();
  }


  template<int dims, int sp_dims>
  inline TensorBlossom<dims, sp_dims> :: TensorBlossom(int dg[dims], int sz[dims]) : domain_computed(false) 
  { 
    copy(dg, dg+dims, deg);
    ctl_msh.Resize(sz);

    unset_periodic();
    setKt();
    set2bsp_srf_type();
    init_iso_generator();
  }


  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: generate_iso() // compute isos in all end points of currrent sbudivided bzr
  {
    iso_generator();
  }
  
    

  template<int dims, int sp_dims>
  template<int from_sp_dims>
  inline TensorBlossom<dims, sp_dims> :: TensorBlossom(int method, SubMFType const*from_mf1, SubMFType const*from_mf2) : domain_computed(false)
  {
    switch(method)
    {
      case revolution: construct_revolution_manifold(*from_mf1); break;
      case ruled: construct_ruled_manifold(*from_mf1, *from_mf2); break;
    }
    init_iso_generator();
  }
  
  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: 
  template<int from_sp_dims>
  construct_ruled_manifold(SubMFType const& from_mf1, SubMFType const & from_mf2)
  {
    SubMFType  mf1(from_mf1), mf2(from_mf2);

    cout << "mf1 is \n" << mf1 << endl;
    cout << "mf2 is \n" << mf2 << endl;
    

    mf1.Convert2OpenEnd(); mf2.Convert2OpenEnd();
    cout << "mf1 is \n" << mf1 << endl;
    cout << "mf2 is \n" << mf2 << endl;
    
    for(int i=0; i<dims-1; i++)
    {
      if(mf1.deg[i] > mf2.deg[i])
        mf2.RaiseDegree(mf1.deg[i] - mf2.deg[i], i);
      else if(mf1.deg[i] < mf2.deg[i])
        mf1.RaiseDegree(mf2.deg[i] - mf1.deg[i], i);
    }
    cout << "\n\nAfter degree raising\n";
    
    cout << "mf1 is \n" << mf1 << endl;
    cout << "mf2 is \n" << mf2 << endl;
    
    for(int i=0; i<dims-1; i++)
    {
      multiset<double> kv1(mf1.ctl_kv[i].begin(), mf1.ctl_kv[i].end()), kv2(mf2.ctl_kv[i].begin(), mf2.ctl_kv[i].end());
      multiset<double> kv12, kv21;    // kv12 = kv1 - kv2, ...
      set_difference( kv1.begin(), kv1.end(), kv2.begin(), kv2.end(), inserter(kv12, kv12.begin()) );
      set_difference( kv2.begin(), kv2.end(), kv1.begin(), kv1.end(), inserter(kv21, kv21.begin()) );
      for(multiset<double>::iterator itr = kv12.begin(); itr != kv12.end(); ++itr)
      {
        mf2.Insert(*itr, i);
      }
      cout << "\n\nAfter refine\n";
      cout << "mf2 is \n" << mf2 << endl;
      for(multiset<double>::iterator itr = kv21.begin(); itr != kv21.end(); ++itr)
      {
        mf1.Insert(*itr, i);
      }
      cout << "now mf1 is \n" << mf1 << endl;
      
      set_union( kv1.begin(), kv1.end(), kv2.begin(), kv2.end(), back_insert_iterator< vector<double> > (ctl_kv[i]) );
    }
    ctl_kv[dims-1].push_back(0); ctl_kv[dims-1].push_back(1);                    // knot vector for the ruled direction.

    int sz[dims]; mf1.ctl_msh.GetSize(sz); sz[dims-1] = 2; ctl_msh.Resize(sz);   // set ctl_msh size.

    slc_iterator itr = ctl_msh.SliceIterate(dims-1);                             // and then copy mf1 and mf2 as its two slices.
    copy(mf1.ctl_msh.Begin(), mf1.ctl_msh.End(), itr.Begin());
    copy(mf2.ctl_msh.Begin(), mf2.ctl_msh.End(), (++itr).Begin());
    
    copy(&mf1.deg[0], &mf1.deg[0] + dims - 1, deg); deg[dims-1] = 1;                     // Set deg in dims-1 to be linear and Copy the rest.
    copy(&mf1.periodic[0], &mf1.periodic[0] + dims - 1, periodic); periodic[dims-1] = 0; // UnSet periodic in dims-1 and Copy the rest.
    rational = from_mf1.rational || from_mf2.rational;                                   // rational if either of from_mf is.

    set2bsp_srf_type();
  }
  
    
  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: 
  template<int from_sp_dims>
  construct_revolution_manifold(TensorBlossom<dims-1, from_sp_dims> const & _rev_from)
  {
    assert(sp_dims==4 && sp_dims>=from_sp_dims && dims>1);

    SubMFType & rev_from = (SubMFType&) _rev_from;
    
    int sz[dims]; rev_from.ControlMesh().GetSize(sz); sz[dims-1] = 8;        // sz[dims-1] is for ctrl point for the circle, the added manifold dimension.
    ctl_msh.Resize(sz);
    copy(rev_from.ctl_msh.Begin(), rev_from.ctl_msh.End(), ctl_msh.LinearIterate(dims-1, 0, 1));
    
    if(from_sp_dims == sp_dims)
    {
      assert(rev_from.IsRational());
    }
    else
    {
      assert(from_sp_dims + 1 == sp_dims || from_sp_dims + 2 == sp_dims);
      
      for(linear_iterator itr = ctl_msh.LinearIterate(dims-1, 0, 1); !itr.PastEnd(); ++itr)
      {
        if(from_sp_dims == 3 && rev_from.IsRational()) // 1D or 2D R in P2.
        {
          (*itr)[sp_dims-1] = (*itr)[2];               // [2] is originally weight component.
          (*itr)[2] = 0;  
        }
        else                                            // 1D or 2D P in E2 or E3 
          (*itr)[sp_dims-1] = 1.0;                      // add weight component.
      }
    }
    /* rotate each original ctl point around Y-axis to generate new contol points. */
    for(row_iterator r_itr = ctl_msh.RowIterate(dims-1); !r_itr.PastEnd(); ++r_itr)
    {
      element_iterator itr0 = r_itr.Begin();
      element_iterator itr = r_itr.Begin();

      point_t& p = *itr0;
      double w = sqrt(2.0) / 2;
      double ang = 45;

      (*++itr) = rY(ang)(p).Syw(w);

      for(++itr, ang += 45; itr != r_itr.End(); ++itr, ang += 45)
      {
        (*itr) = rY(ang)(p);
        ang += 45;
        (*++itr) = rY(ang)(p).Syw(w);
      }
    }
    for(int i=0; i<dims-1; i++)                                    // copy the original knot vectors.
    {
      copy(rev_from.ctl_kv[i].begin(), rev_from.ctl_kv[i].end(), back_insert_iterator< vector<double> >(ctl_kv[i]) );
    }
    for(int i=0; i<4; i++)                                         // set knot vector along the parallel to be 0,0,1,1,2,2,3,3,4.
    {
      ctl_kv[dims-1].push_back(i); ctl_kv[dims-1].push_back(i);
    }
    ctl_kv[dims-1].push_back(4);

    copy(rev_from.deg, rev_from.deg + dims-1, deg);                 // copy degrees.
    copy(rev_from.periodic, rev_from.periodic + dims-1, periodic);  // copy periodic condition.

    deg[dims-1] = 2;                                                // Set deg in dims-1 to be quadratic and,
    periodic[dims-1] = true;                                        // of course periodic and
    IsRational() = true;                                            // rational.

    set2bsp_srf_type();
  }
  
      
    
  template<int dims, int sp_dims>
  inline TensorBlossom<dims, sp_dims>& TensorBlossom<dims, sp_dims>::Transform(xchen::Transform const& trans)
  {
    ctl_msh.Transform(trans);
    clean_sbd();
    return *this;
  }
  
    
  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: SetUniformFloatingEndKnots()
  {
    for(int i=0; i<dims; i++)
    {
      generate( ctl_kv[i].begin(), ctl_kv[i].end(), iota<double>(0.0) );
    }
  }
  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: SetUniformOpenEndKnots()
  {
    for(int i=0; i<dims; i++)
    {
      vector<double>::iterator itr = ctl_kv[i].begin() + deg[i];

      fill(ctl_kv[i].begin(), itr, 0.0);

      double u = 1.0;
      for(; itr != ctl_kv[i].end() - deg[i]; ++itr, ++u)
        *itr = u;
      
      fill(itr, ctl_kv[i].end(), u);
    }
  }

  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: init_sbd()
  {
    if(sbd_initialized) return;
    
    BaseType::init_sbd(); 
    if(!domain_computed) computeSegmentDomains(); 

    for(int i=0; i<dims; i++)
    {
      init_sbd_knotvector(i);
      convert2Bzr(i, sbddata);
    }
    compact_sbd_kv();
    sbd_initialized = true;
  }


  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: compact_sbd_kv()
  {
    for(int i=0; i<dims; i++)
    {
      vector<double>::iterator in_itr = sbd_kv[i].begin() + deg[i]-1, out_itr = sbd_kv[i].begin();
      for(; in_itr != sbd_kv[i].end()+deg[i]-1; in_itr += deg[i], ++out_itr)
      {
        *out_itr = *in_itr;
      }
      sbd_kv[i].erase(out_itr, sbd_kv[i].end());
    }
  }
  
  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: subdivide_kv()
  {
    for(int i=0; i<dims; i++)
    {
      vector<double>& kv = (vector<double>&) sbd_kv[i];
      vector<double> old_kt(kv);
      kv.resize(2*old_kt.size()-1);
      vector<double>::iterator itr1 = old_kt.begin(), itr2 = itr1 + 1, itr = kv.begin()-1;
      for(; itr2 != old_kt.end(); ++itr1, ++itr2)
      {
        *(++itr) = *itr1; *(++itr) = (*itr1 + *itr2)/2;
      }
      *(++itr) = *itr1; 
    }
  }
  

  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: init_sbd_knotvector(int i)
  {
    sbd_kv[i] = ctl_kv[i];
    if(periodic[i]) 
      periodic2floating(i, sbddata);
  }

  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: periodic2floating(int i, int applyon)
  {
    knot_t& kv = cur_kv(applyon);
    mesh_t& msh = cur_msh(applyon); 

    double T = kv[i].back() - kv[i].front();        // modula of T
    vector<double> newkv;
    back_insert_iterator< vector<double> > bitr(newkv);
        
    transform(kv[i].end()-deg[i], kv[i].end()-1, bitr, bind2nd(minus<double>(), T));     // copy the last deg-1 upto back() ele(-T).
    copy(kv[i].begin(), kv[i].end(), bitr);                                              // copy the whole knots
    transform(kv[i].begin()+1, kv[i].begin()+deg[i], bitr, bind2nd(plus<double>(), T));  // copy the first deg-1 ele(+T).

    kv[i] = newkv;

    int extendRL[2*dims];                                                // will append front deg-1 and back 1 extra pnts.
    bzero(extendRL, sizeof(int) * 2 * dims);
    extendRL[2*i] = 1; extendRL[2*i+1] = deg[i] - 1;
    msh.Extend(extendRL);
    
    array_copy(i, deg[i]-1,  msh, msh.GetSize(i) - deg[i],  msh, 0);     // fill the left expanded deg-1 elements from the right eles(before expanded).
    array_copy(i, 1,  msh, deg[i]-1,  msh, msh.GetSize(i) - 1);          // fill the right expanded 1 element from the first(before expanded) 

    if(applyon==ctldata)
      periodic[i] = 0;
  }
    

  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: identifyEndPntsForPeriodicEnd(int i, int applyon)
  {
    if(periodic[i])
    {
      mesh_t msh = cur_msh(applyon);
      if(msh.GetSize(i) > 1)
        for(row_iterator itr = msh.RowIterate(i); !itr.PastEnd(); ++itr)
          *(itr.End()-1) = *(itr.Begin());
    }
  }
  
  
  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: convert2Bzr(int i, int applyon)
  {
    convert2OpenEnd(i, applyon);

    if(dmn[i].size() > 2)
    {
      vector<double>::iterator itr = dmn[i].begin()+1, end_itr = dmn[i].end() - 1;
      for(; itr != end_itr; ++itr)
      {
        while(insert(*itr, i, applyon));
      }
    }
  }

  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: convert2OpenEnd(int i, int applyon)
  {
    knot_t& kv = cur_kv(applyon);
    mesh_t& msh = cur_msh(applyon); 

    if(!domain_computed) computeSegmentDomains(); 
    
    if(applyon==ctldata && periodic[i])
    {
      periodic2floating(i, ctldata);
    }
    while( insert(dmn[i].front(), i, applyon) ) { popFront(i, applyon);  }
    while( kv[i].front() != dmn[i].front() )     popFront(i, applyon);

    while( insert(dmn[i].back(), i, applyon) )  { popBack(i, applyon);}
    while( kv[i].back() != dmn[i].back() )       popBack(i, applyon);

    identifyEndPntsForPeriodicEnd(i, applyon);
  }
  
  
  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: computeSegmentDomains() 
  { 
    for(int i=0; i<dims; i++)
    {
      if(!dmn[i].size())
      {
        vector<double>::iterator itr1 = periodic[i]? ctl_kv[i].begin() : ctl_kv[i].begin() + deg[i] - 1;
        dmn[i].push_back(*itr1);
      
        vector<double>::iterator itr2 = itr1 + 1;
        for(; itr2 != (periodic[i]? ctl_kv[i].end() : ctl_kv[i].end() - deg[i] + 1); ++itr1, ++itr2)
        {
          if(*itr1 != *itr2)
          {
            dmn[i].push_back(*itr2);
          }
        }
      }
    }
    domain_computed = true;
  }
  
  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: RaiseDegree(int inc, int dir) 
  {
    Convert2Bzr(dir);

    mesh_t old_msh(ctl_msh);
    ctl_msh.ExtendAtDirection(dir, 0, inc * (ctl_msh.GetSize(dir) - 1) / deg[dir]);
    fill(ctl_msh.Begin(), ctl_msh.End(), point_t(0.0));

    int n = deg[dir],  new_deg = n+inc;
    
    linear_iterator* out_itr = new linear_iterator[new_deg];                // iterate each pnt in a bzr polygon for the new mesh.

    slc_iterator itr = old_msh.SliceIterate(dir); 
    for(int poly = 0; poly < (old_msh.GetSize(dir) - 1) / n; poly++)        // iterate all bzr control 1D polygon along i direction.
    {
      for(int i=0; i < new_deg; i++)
        out_itr[i] = ctl_msh.LinearIterate(dir, poly * new_deg + i, 1);     // used to iterate all slices other than the right end one.
      
      for(int i=0; i <= n; ++i, ++itr)                                      // iterate all pnts within the polygon.
      {
        for(int j=0; j <= inc; j++)                                         // b<0^(n-i)1^i> contributing to b<0^(n-i + inc-j)1^(i + j)>
        {
          if( (!i && !j) || (i==n && j==inc) ) continue;                    // can only contribute to the two end pnts.
          
          transform(itr.Begin(), itr.End(), out_itr[i+j].Begin(), out_itr[i+j].Begin(), 
                    Scale1stAdd<point_t>( C(j, i+j) * C(inc-j,  n-i + inc-j) ));
        }
      }
      --itr;                                                                // the last pnt of this poly is also the first pnt of the next poly.
    }
    delete [] out_itr;
    
    transform( ctl_msh.Begin(), ctl_msh.End(), ctl_msh.Begin(), ScaleAll<point_t>( 1.0/C(inc, n + inc) ) );   // make it affine combination.
    slc_iterator itr1 = old_msh.SliceIterate(dir), itr2 = ctl_msh.SliceIterate(dir);                  // and copy all the end points of bzr polys.
    for(; !itr1.PastEnd(); itr1 += n, itr2 += new_deg)
    {
      copy(itr1.Begin(), itr1.End(), itr2.Begin());
    }

    deg[dir] += inc;
    for(vector<double>::iterator itr = dmn[dir].begin(); itr != dmn[dir].end(); ++itr)
      for(int i=0; i<inc; i++)
        ctl_kv.Insert(*itr, dir);
    
    clean_sbd();
    return;
  }

  template<int dims, int sp_dims>
  inline bool TensorBlossom<dims, sp_dims> :: is_valid_range(IdxRange const& range, int applyon)
  {
    return range.s+1 < range.e && range.s >=0 && range.e <= cur_msh(applyon).GetSize(range.d);
  }
  
  template<int dims, int sp_dims>
  inline bool TensorBlossom<dims, sp_dims> :: Insert(double u,int i)
  { 
    return !periodic[i] && insert(u,i,ctldata); 
  }
  template<int dims, int sp_dims>
  bool TensorBlossom<dims, sp_dims> :: insert(double u, int dir, int applyon)
  {
    assert(deg[dir]);    

    knot_t& kv = cur_kv(applyon);
    mesh_t& msh = cur_msh(applyon);
    
    TensorBlossom<dims, sp_dims>*This = (TensorBlossom<dims, sp_dims>*)this;
    
    int repeat;                                          // The multiplicity of the u value in current knot vector.
    int idx = kv.UpperBound(u, repeat, dir) - deg[dir];  // First idx such that kv[dir][idx] > u - deg = 1st pnt of the seg where u is defined.
    IdxRange range(dir, idx, idx + deg[dir] + 1 - repeat);

    if( ! is_valid_range(range, applyon) ) return false;

    vector<point_t> first_pnts;      // save these first pnts as they are going to be overwritten by insert().
    copy( msh.LinearIterate(dir,idx,1).Begin(), msh.LinearIterate(dir,idx,1).End(), back_insert_iterator<vector<point_t> >(first_pnts) );

    This->blend(range, u, applyon);           // blending ctl points. the first pnt is overwritten, and the last pnt is kept.
    kv.Insert(u, dir);               // insert the u knot value.

    if(dims==1)               
      msh.Insert(msh.Begin() + idx, first_pnts.front());
    else 
    {
      msh.Insert(dir, idx);
      copy(first_pnts.begin(), first_pnts.end(), msh.LinearIterate(dir, idx, 1));   // now restore the first pnt.
    }
    return true;
  }
  

  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: DecVarity(double u, int dir)
  {
    for(row_iterator itr = ctl_msh.RowIterate(dir); !itr.PastEnd(); ++itr)
    {
      insert(dir, itr.Begin(), itr.End(), itr.Begin(), 0, u);
    }

    deg[dir]--;
    ctl_msh.ResizeOneDirection(dir, ctl_msh.GetSize(dir)-1);
    ctl_kv[dir].pop_back();
    ctl_kv[dir].erase(ctl_kv[dir].begin());
  }


  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: blend(IdxRange const& range, double u, int applyon) 
  {
    knot_t& kv = cur_kv(applyon);
    mesh_t& msh = cur_msh(applyon);

    int d = range.d, s = range.s, e = range.e;
    int i = s,  j = s + deg[d];
    slc_iterator 
      input1 = msh.SliceIterate(d, s, e-1), 
      input2 = msh.SliceIterate(d, s+1, e), 
      output = msh.SliceIterate(d, s, e-1); 

    for(; !input1.PastEnd(); ++input1, ++input2, ++output, ++i, ++j)
    {
      transform(input1.Begin(), input1.End(), input2.Begin(), output.Begin(), 
                AffineCombine<point_t>( (u - kv[d][i]) / (kv[d][j] - kv[d][i]) ));
    }
  }
  
  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: Subdivide(int depth) const
  {
    TensorBlossom<dims, sp_dims>*This = ((TensorBlossom<dims, sp_dims>*)this);
    This->init_sbd();
    mesh_t& msh = This->sbd_msh;

    msh.SetN(depth);
    while(depth--)
    {
      msh.InsertMiddlePoints(false);
      for(int i=0; i<dims; i++)
      {
        for(slc_iterator itr = msh.SliceIterate(i); !(itr+1).PastEnd(); itr += 2*deg[i])           // iterate all bzr control 1D polygon along i direction.
        {
          slc_iterator left_itr0(itr), middle_itr0(itr + 1), right_itr0(itr + 2); 

          for(int j=0; j < deg[i]; j++, ++left_itr0, ++right_itr0, ++middle_itr0)                  // de-castejau algo from level 0 to level deg[i]-1.
          {
            slc_iterator left_itr(left_itr0), right_itr(right_itr0), middle_itr(middle_itr0);

            for(int k=0; k<deg[i]-j;  k++, left_itr += 2, right_itr += 2, middle_itr += 2)         // de-castejau algo level j.
            {
              transform(left_itr.Begin(), left_itr.End(), right_itr.Begin(), middle_itr.Begin(), Average<point_t>() );
            }
          }
        }
      }
      This->subdivide_kv();
    }
    This->compute_1st_order_derivative();
    //    This->generate_iso();
  }

  template<int dims, int sp_dims>
  void TensorBlossom<dims, sp_dims> :: compute_1st_order_derivative()
  {
    if(dims>2) return;
    
    int sz[dims];
    sbd_msh.GetSize(sz);
    for(int i=0; i<dims; i++) sz[i] = (sz[i]-1) / deg[i] + 1; 
    for(int i=0; i<dims; i++) Dmsh[i].Resize(sz);

    for(int i=0; i<dims; i++)
    {
      row_iterator msh_row_itr = sbd_msh.RowIterate(i), Dmsh_row_itr = Dmsh[i].RowIterate(i);

      for(; !Dmsh_row_itr.PastEnd(); ++Dmsh_row_itr, msh_row_itr += deg[(i+1)%dims])   // only for dims == 2;  need more general algo here.
      {
        element_iterator msh_ele_itr1 = msh_row_itr.Begin(), msh_ele_itr2 = msh_row_itr.Begin() + 1, Dmsh_ele_itr = Dmsh_row_itr.Begin();

        for(; msh_ele_itr2 < msh_row_itr.End(); msh_ele_itr1 += deg[i], msh_ele_itr2 += deg[i], ++Dmsh_ele_itr)
        {
          *Dmsh_ele_itr = *msh_ele_itr2 - *msh_ele_itr1;       // need to be scaled by domain interval and deg.
        }
        *Dmsh_ele_itr = *msh_ele_itr1 - *(msh_ele_itr1-1);
        if(dims==1) break;
      }
    }
    if(IsRational())
    {
      for(int i=0; i<dims; i++)      // compute 1st order partial derivative along direction i.
      {
        linear_iterator Dmsh_itr = Dmsh[i].Begin();
        for(row_iterator msh_row_itr = sbd_msh.RowIterate(0); !msh_row_itr.PastEnd(); msh_row_itr += deg[1]) //assume dims==2 now.
        {
          int count = 0;
          for(element_iterator msh_itr = msh_row_itr.Begin(); msh_itr <  msh_row_itr.End(); msh_itr += deg[0], ++Dmsh_itr)
          {
            double w = (*msh_itr)[3], w1 = (*Dmsh_itr)[3];          // last components of the point and the derivative
            Vector<double, sp_dims-1> p(*msh_itr), p1(*Dmsh_itr);   // the point and the derivative without last component.
            *Dmsh_itr = (w * p1 - w1 * p) / (w * w);
          }
          if(dims==1) break;
        }
      }
    }
    if(dims==2)
    {
      Nmsh.Resize(sz);
      transform(Dmsh[1].Begin(), Dmsh[1].End(), Dmsh[0].Begin(), Nmsh.Begin(), cross_product());
    }
  }
  
  template<int dims, int sp_dims>
  inline void TensorBlossom<dims, sp_dims> :: draw_iso()
  {
    draw_iso_functor<dims, sp_dims>(*this)();
  }
  
        
  template<int dims, int sp_dims>
  ostream& operator<< (ostream& os, TensorBlossom<dims, sp_dims> const& _blsm)
  {
    TensorBlossom<dims, sp_dims>& blsm = (TensorBlossom<dims, sp_dims> &) _blsm;
    
    os << "The " << (blsm.IsRational()? "rational" : "polynomial") << " blossom has,\nControl Mesh = \n" << blsm.ctl_msh << endl;
    os << "Degrees = \t[";
    for(int i=dims-1; i; i--)
      os << blsm.deg[i] << " ";
    os << blsm.deg[0] << "]\n\n";
    os << "Knot vectors =\t";
    for(int i=dims-1; i>=0; i--) 
    {
      os << "[" << (blsm.IsPeriodicEnd(i)? "Periodic: " : ""); 
      copy(blsm.ctl_kv[i].begin(), blsm.ctl_kv[i].end()-1, ostream_iterator<double>(os, " ")); os << blsm.ctl_kv[i].back() << "]\t";
    }
    os << endl;

    return os;
  }




  template<int dims, int sp_dims>
  inline TensorBlossom<dims, sp_dims> operator*(Transform const& trans, TensorBlossom<dims, sp_dims>const& ar)
  {
    TensorBlossom<dims, sp_dims> ret(ar);
    ret.Transform(trans);
    return ret;
  }
  
    
}//end namespace xchen


#endif

---