---

SubdivisionManifold.H

Go to the documentation of this file.

#ifndef _SUBDIVISIONMANIFOLD_H
#define _SUBDIVISIONMANIFOLD_H


#include <TensorManifold.H>
#include <arglst.H>
#include <GL/gl.h>
#include <xstl.H>


namespace xchen
{
  template<int dims = 2, int sp_dims = 3>
  class SubdivisionManifold : public TensorManifold<dims, sp_dims>
  {
      int len[dims];
      int left_ec[dims], right_ec[dims];             // floating, open, or periodic end conditions.
      enum {open_end, floating_end, periodic_end};

      def_type_name_for_tensor_manifold(dims, sp_dims);
      
  public:
      SubdivisionManifold(int deg1, ...);            
      SubdivisionManifold(int* degsz);               

      void Subdivide(int depth=1)const;              

      void SetOpenCondition(int* l_ec, int* r_ec);   
      void SetOpenCondition(int ec);                 
       
  private:
      void set_default_end_condition()               { for(int i=0; i<dims; i++) left_ec[i] = right_ec[i] = open_end; }   
      void convert2floatingends();              
      void extrapolate2floating(int dir, bool left);
      void init_subdivision_data()                   { init_sbd_mesh();   convert2floatingends(); }

  public:
      void bzr2floatingEndBspline(int dir);
  };



  template<int dims, int sp_dims>
  inline SubdivisionManifold<dims, sp_dims> :: SubdivisionManifold(int arg1, ...) 
  { 
    int args[2*dims];
    get_reverse_ints_va(args, 2*dims); 
    for(int i=dims-1; i>=0; i--)
    {
      deg[i] = args[2*i+1];
      len[i] = args[2*i];
    }
    ctl_msh.Resize(len);
    set_default_end_condition();
  }
  template<int dims, int sp_dims>
  inline SubdivisionManifold<dims, sp_dims> :: SubdivisionManifold(int* degsz)
  { 
    for(int i=0; i<dims; i++)
    {
      deg[i] = degsz[2*i];
      len[i] = degsz[2*i+1];
    }
    msh.Resize(len);
    set_default_end_condition();
  }


  template<int dims, int sp_dims>
  void SubdivisionManifold<dims, sp_dims> :: SetOpenCondition(int* l_ec, int* r_ec)  
  {
    memcpy(left_ec, l_ec, sizeof(int) * dims);
    memcpy(right_ec, r_ec, sizeof(int) * dims);
  }
  template<int dims, int sp_dims>
  void SubdivisionManifold<dims, sp_dims> :: SetOpenCondition(int ec)
  {
    fill(left_ec, left_ec+dims, ec);
    fill(right_ec, right_ec+dims, ec);
  }
  template<int dims, int sp_dims>
  void SubdivisionManifold<dims, sp_dims> :: convert2floatingends()
  {
    for(int i=0; i<dims; i++)
    {
      if(left_ec[i] == open_end && right_ec[i] == open_end && msh.GetSize(i) == deg[i] + 1)
      {
        cerr << "Subdivide() can only apply to tensor mesh of same degree in all directions.\n"; 
        bzr2floatingEndBspline(i);
      }
      
      else 
      {
        if(left_ec[i] == open_end)
          extrapolate2floating(i, true);
        if(right_ec[i] == open_end)
          extrapolate2floating(i, false);
      }
    }
  }
  

  template<typename T>
  struct affine_combine_inc_t1 : public binary_function<T, T const&, T const&>
  {
      affine_combine_inc_t1(double t_1, double t_2) : t1(t_1), t2(t_2) { }
      T operator()(T const& p1, T const& p2) { double t11 = t1++; return (p1 * t11 + p2 * t2) / (t11 + t2); }
      double t1, t2;
  };
  template<int dims, int sp_dims>
  void SubdivisionManifold<dims, sp_dims> :: extrapolate2floating(int dir, bool left)
  {

    int d = deg[dir];
    for(row_iterator rw_itr = msh.RowIterate(dir); !rw_itr.PastEnd(); ++rw_itr)
    {
      int t10 = 2, t2 = -1;
      for(int i=1; i < d; i++, t10++, t2--)
      {
        int t1 = t10;
        if(left) 
          transform(rw_itr.Begin(), rw_itr.Begin()+(d-i), rw_itr.Begin()+1, rw_itr.Begin(), affine_combine_inc_t1<point_t>(t1, t2));
        else
          transform(rw_itr.Rbegin(), rw_itr.Rbegin()+(d-i), rw_itr.Rbegin()+1, rw_itr.Rbegin(), affine_combine_inc_t1<point_t>(t1, t2));
      }
    }
  }
  
  template<int dims, int sp_dims>
  void SubdivisionManifold<dims, sp_dims> :: bzr2floatingEndBspline(int dir)
  {
    cout << "\nctl mesh is\n" << msh << endl;
    int n = deg[dir];
    point_t deg2_mid_pnt;
    for(row_iterator rw_itr = msh.RowIterate(dir); !rw_itr.PastEnd(); ++rw_itr)
    {
      if(n==2)
        deg2_mid_pnt = *(rw_itr.Begin() + 1);

      double *kv = new double[2 * n];
      fill(kv, kv+n, 0); fill(kv+n,kv+2*n, 1);
      TensorArray<point_t, 1> new_pnts(n+1);
      {
        Blossom<sp_dims> blsm(n, n);
        copy(rw_itr.Begin(), rw_itr.End()-1, blsm().Begin());
        blsm.KnotVector() = kv;
        for(int i = 0; i < n-1; i++)               // evaluate n-1 times with u = -1, -2, ..., -(n-1).
        {
          blsm.Evaluate( -(i+1) );
          if( i > n-4 )                            // last two evaluations of i=n-3, i=n-2. Total evaluation times are n-2, and n-1.
            new_pnts[n-2-i] = blsm()[n-2-i];
          else                                     // have evaluated(to left) at u=-1,.. with i+1<= n-3 times so far, 
          {
            Blossom<sp_dims> blsm2(blsm);
            blsm2.PopFont();
            for(int j=0; j<n-i-3; j++)             // need to evaluate(to right) n-i-3 times, with u=2, 3,..
              blsm2.Evaluate(2+j);
            new_pnts[n-i-2] = blsm2()[0];          // blossom value at knot values starting at -(i+1), which is -(i+1)-(1-n)'th pnt.
          }
        }
      }
      {
        Blossom<sp_dims> blsm(n, n);
        copy(rw_itr.Begin()+1, rw_itr.End(), blsm().Begin());
        blsm.KnotVector() = kv+1;
        for(int i = 0; i < n-1; i++)               // evaluate n-1 times with u = 2, 3, ..., n.
        {
          blsm.Evaluate(2 + i);
          if( i > n-4 )                            // last two evaluations of i=n-3, i=n-2.
            new_pnts[2+i] = blsm()[0];
        }
      }
      copy(new_pnts.Begin(), new_pnts.End(), rw_itr.Begin());
      if(n==2)
        *(rw_itr.Begin() + 1) = deg2_mid_pnt;

      cout << "\nctl mesh is\n" << msh << endl;
    }
  }
  
  template<int dims, int sp_dims>
  void SubdivisionManifold<dims, sp_dims> :: Subdivide(int depth) const
  {
    int d = deg[0];
    for(int i=1; i<dims; i++)
      if( d != deg[i] )
      {
        cerr << "Subdivide() can only apply to tensor mesh of same degree in all directions.\n", exit(0);
      }
    
    init_sbd_msh();
    
    if(sbdl_msh.Empty())
    { 
      original_msh = msh;
      convert2floatingends();
    }

    msh.SetN(depth);
    while(depth--)
    {
      msh.InsertMiddlePoints();
      for(int i=2; i <= deg[0]; i++)
        msh.Averaging();
    }
  }

  

  
  
}//end namespace xchen

#endif

---