include/shape/bcell3d_algebraic_surface.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a  g i x
 *****************************************************************************
 * bcell of algebraic surface
 * 2008-03-20
 * Bernard Mourrain & Julien Wintz
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/

# ifndef shape_bcell3d_algebraic_surface_hpp
# define shape_bcell3d_algebraic_surface_hpp

# include <realroot/Interval.hpp>
# include <realroot/ring_bernstein_tensor.hpp>
# include <realroot/sign_variation.hpp>
# include <realroot/Seq.hpp>
# include <shape/bounding_box.hpp>
# include <shape/algebraic_surface.hpp>
# include <shape/bcell3d.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/marching_cube.hpp>

# define TMPL template<class C, class V>
# define TMPL1 template<class C>
# define SELF bcell3d_algebraic_surface<C,V>
# undef Topology
//====================================================================
namespace mmx { namespace shape {
//====================================================================

// TMPL struct topology;
// TMPL struct tpl3d;
// TMPL struct bcell3d_algebraic_surface;

TMPL class topology;
TMPL class tpl3d;
TMPL class bcell3d_algebraic_surface;

//--------------------------------------------------------------------
struct bcell3d_algebraic_surface_def {}; 

template<> struct use<bcell3d_algebraic_surface_def>
//   :public use<topology_def> 
//   ,public use<algebraic_surface_def> 
{
  typedef Interval< double >                    Scalar;
  typedef polynomial< double, with<Bernstein> > Polynomial;
  //  typedef use<bcell3d_def>::Cell3d           Cell3d;
  //  typedef use<bcell3d_def>::BoundingBox      BoundingBox;
  //  typedef bcell3d_algebraic_surface<double,double> Cell3dAlgebraicSurface;
}; 

//--------------------------------------------------------------------
template<class C, class V=default_env>
class bcell3d_algebraic_surface : public bcell3d<C,V> {
public:

  typedef REF_OF(V)  Ref;
  typedef tpl3d<C,V>                          Topology3d;
  typedef topology<C,REF_OF(V)>               Topology;

  typedef typename Topology::Point         Point;  
  typedef typename Topology::Edge          Edge;
  typedef typename Topology::Face          Face;

  typedef bcell3d<C,V>                      Cell;
  typedef cell3d<C,V>                       CellBase;
  typedef bcell3d<C,V>                      Cell3d;

  typedef bounding_box<C,Ref>              BoundingBox;
  
  typedef algebraic_surface<C,Ref>         AlgebraicSurface;
  
  typedef typename use<bcell3d_algebraic_surface_def,C,V>::Polynomial  Polynomial;

  //  typedef polynomial< Interval<C>, with<Bernstein> > Polynomial;
  //  typedef polynomial< Scalar, with<Bernstein> > Polynomial;
  
  typedef solver_implicit<C,V>             Solver;
  
  bcell3d_algebraic_surface(const SELF&);
  bcell3d_algebraic_surface(const Polynomial&, const BoundingBox& );
  bcell3d_algebraic_surface(char*, const BoundingBox& );
  
  bcell3d_algebraic_surface(const AlgebraicSurface&, const BoundingBox&);
  bcell3d_algebraic_surface(AlgebraicSurface *, const BoundingBox&);

  virtual bool is_regular(void);
  virtual bool is_active (void) const;

  virtual bool is_intersected (void) {return true;}

  virtual bool insert_regular (Topology*) ;
  virtual bool insert_singular(Topology*) ;

  virtual void polygonise (Topology3d*) ;

  virtual bool topology_regular(Topology*) ;
  
  template<class CELL> 
          void split    (CELL*& left, CELL*& right, int v, double s);
  virtual void subdivide(CellBase*& left, CellBase*& right, int v, double s);
  
  int          mc_index(void)      const ;
  bool         edge_point(Point**, int ) const ;

  Point        center() const { return Point((this->xmin()+this->xmax())/2, (this->ymin()+this->ymax())/2,
                                       (this->zmin()+this->zmax())/2); }

  C            center_value() { return m_polynomial(0.5,0.5,0.5); }

  void         insert(Point* p) { m_points<< p; }
  void         insert(Face*  p) { m_faces<< p;  }

  double       vertex_eval(unsigned sx, unsigned sy, unsigned sz) const;
  const Polynomial&   equation() const {return m_polynomial;}
  Polynomial   get_polynomial() const {return m_polynomial;}
  
public:
  Polynomial        m_polynomial;
  Point*            m_center;
  C                 m_center_value;
  int               m_idx;

  Seq<Point*>       m_points;
  Seq<Face*>        m_faces;

private:
  void resolve_conflict(Seq<Cell3d*> &neighb, double side, int cc);
};

//====================================================================
TMPL
SELF::bcell3d_algebraic_surface(const SELF& s)
  : Cell3d((BoundingBox)s), m_polynomial(s.equation()), m_center(0)
{
  
}
//--------------------------------------------------------------------
TMPL
SELF::bcell3d_algebraic_surface(const Polynomial& pol, const BoundingBox& bx)
  : Cell3d(bx), m_polynomial(pol), m_center(0)
{

}
//--------------------------------------------------------------------
TMPL
SELF::bcell3d_algebraic_surface(char* pol, const BoundingBox& bx)
  : Cell3d(bx), m_polynomial(pol), m_center(0)
{
}
//--------------------------------------------------------------------
TMPL
SELF::bcell3d_algebraic_surface(const AlgebraicSurface& sf, const BoundingBox& b): Cell3d(b), m_center(0)
{
  Seq<mmx::GMP::rational> bx;
  bx<<as<mmx::GMP::rational>(b.xmin());
  bx<<as<mmx::GMP::rational>(b.xmax());
  bx<<as<mmx::GMP::rational>(b.ymin());
  bx<<as<mmx::GMP::rational>(b.ymax());
  bx<<as<mmx::GMP::rational>(b.zmin());
  bx<<as<mmx::GMP::rational>(b.zmax());
  
  Polynomial tmp;
  tensor::bernstein<mmx::GMP::rational> polq(sf.equation().rep(),bx);
  let::assign(tmp.rep(),polq);
  m_polynomial=tmp;

//     BoundingBox* bbx = sf->boundingBox();
//     Polynomial ft, bk;
//     tensor::face(ft,m_polynomial, 2, 1);
//     solver_implicit::edge_point(n_intersections, ft, north_edge, bbx);
//     solver_implicit::edge_point(s_intersections, ft, south_edge, bbx);
//     solver_implicit::edge_point(w_intersections, ft, west_edge , bbx);
//     solver_implicit::edge_point(e_intersections, ft, east_edge , bbx);

//     tensor::face(bk, m_polynomial, 2, 0);
//     solver_implicit::edge_point(n_intersections, bk, north_edge, bbx);
//     solver_implicit::edge_point(s_intersections, bk, south_edge, bbx);
//     solver_implicit::edge_point(w_intersections, bk, west_edge , bbx);
//     solver_implicit::edge_point(e_intersections, bk, east_edge , bbx);

//     tensor::face(ft, m_polynomial, 1, 1);
//     solver_implicit::edge_point(w_intersections, ft, west_edge, bbx);
//     solver_implicit::edge_point(e_intersections, ft, east_edge, bbx);

//     tensor::face(bk, m_polynomial, 1, 0);
//     solver_implicit::edge_point(w_intersections, bk, west_edge, bbx);
//     solver_implicit::edge_point(e_intersections, bk, east_edge, bbx);

//     foreach(Point* p, n_intersections) std::cout<<"n "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     foreach(Point* p, s_intersections) std::cout<<"s "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     foreach(Point* p, w_intersections) std::cout<<"w "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     foreach(Point* p, e_intersections) std::cout<<"e "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;

}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_active()  const {
  //BoundingBox* bbx=(BoundingBox*)this;
  
  //    cout<<(BoundingBox*)this<<endl;
  //    cout<<m_polynomial.rep()<<endl;
  //    cout<<has_sign_variation(m_polynomial.begin(),m_polynomial.end())<<endl;
  
  //  Seq<Point*> lp;
  //  Polynomial ft, bk;
  //  tensor::face(ft,m_polynomial, 2, 1);

//     solver_implicit::edge_point(lp, m_polynomial, north_front_edge, bbx);
//     solver_implicit::edge_point(lp, m_polynomial, south_front_edge, bbx);
//     solver_implicit::edge_point(lp, m_polynomial, west_front_edge , bbx);
//     solver_implicit::edge_point(lp, m_polynomial, east_front_edge , bbx);

//     tensor::face(bk, m_polynomial, 2, 0);
//     solver_implicit::edge_point(lp, bk, north_back_edge, bbx);
//     solver_implicit::edge_point(lp, bk, south_back_edge, bbx);
//     solver_implicit::edge_point(lp, bk, west_back_edge , bbx);
//     solver_implicit::edge_point(lp, bk, east_back_edge , bbx);

//     tensor::face(m_polynomial, m_polynomial, 1, 1);
//     solver_implicit::edge_point(lp, ft, west_edge, bbx);
//     solver_implicit::edge_point(lp, ft, east_edge, bbx);

//     tensor::face(bk, m_polynomial, 1, 0);
//     solver_implicit::edge_point(lp, bk, west_edge, bbx);
//     solver_implicit::edge_point(lp, bk, east_edge, bbx);

//    foreach(Point* p, lp) std::cout<<" "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;

  return has_sign_variation(m_polynomial.begin(),m_polynomial.end());
}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_regular() { 
  Polynomial dxf, dyf, dzf;
  Polynomial fp0, fp1;
  Polynomial bk, ft;


  tensor::face(bk, m_polynomial, 2, 0);
  if(Solver::edge_sign_var(bk, Solver::north_back_edge) >1) return false; 
  if(Solver::edge_sign_var(bk, Solver::south_back_edge) >1) return false; 
  if(Solver::edge_sign_var(bk, Solver::east_back_edge ) >1) return false; 
  if(Solver::edge_sign_var(bk, Solver::west_back_edge ) >1) return false; 


  tensor::face(ft, m_polynomial, 2, 1);
  if(Solver::edge_sign_var(ft, Solver::north_front_edge) >1) return false; 
  if(Solver::edge_sign_var(ft, Solver::south_front_edge) >1) return false; 
  if(Solver::edge_sign_var(ft, Solver::east_front_edge ) >1) return false; 
  if(Solver::edge_sign_var(ft, Solver::west_front_edge ) >1) return false; 

  tensor::face(ft, equation(), 1, 1);
  if(Solver::edge_sign_var(ft, Solver::north_west_edge) >1) return false; 
  if(Solver::edge_sign_var(ft, Solver::north_east_edge) >1) return false; 

  tensor::face(bk, equation(), 1, 0);
  if(Solver::edge_sign_var(bk, Solver::south_west_edge) >1) return false; 
  if(Solver::edge_sign_var(bk, Solver::south_east_edge) >1) return false; 

  tensor::diff(dxf.rep(),m_polynomial.rep(),0);
  tensor::diff(dyf.rep(),m_polynomial.rep(),1);
  tensor::diff(dzf.rep(),m_polynomial.rep(),2);

  if(!has_sign_variation(dxf.begin(),dxf.end())) {

    this->m_idx=0;
    return true;
    
    tensor::face(fp0, dyf, 0, 0);
    tensor::face(fp1, dyf, 0, 1);
    if(!has_sign_variation(fp0.begin(),fp0.end()) &&
       !has_sign_variation(fp1.begin(),fp1.end()) ) return true;

    tensor::face(fp0, dzf, 0, 0);
    tensor::face(fp1, dzf, 0, 1);
    if(!has_sign_variation(fp0.begin(),fp0.end()) &&
       !has_sign_variation(fp1.begin(),fp1.end())) return true;

    return false;
  }

  //tensor::diff(dyf.rep(),m_polynomial.rep(),1);
  if(!has_sign_variation(dyf.begin(),dyf.end())) {
    this->m_idx=1;
    return true;

//     tensor::face(fp0, dzf, 0, 0);
//     tensor::face(fp1, dzf, 0, 1);
//     if(!has_sign_variation(fp0.begin(),fp0.end()) &&
//        !has_sign_variation(fp1.begin(),fp1.end())) return true;

//     return false;
  }

  //tensor::diff(dzf.rep(),m_polynomial.rep(),2);
  if(!has_sign_variation(dzf.begin(),dzf.end())) {
    this->m_idx=2;
    return true;
  }

  return false;
}
//--------------------------------------------------------------------
TMPL template<class CELL> void 
SELF::split(CELL*& left, CELL*& right, int v, double c) {
  
  //  cout<<"Split begin"<<endl;
  if(v==0) {
    left = new CELL(m_polynomial, BoundingBox(this->xmin(),c, this->ymin(), this->ymax(), this->zmin(), this->zmax()));
    right= new CELL(m_polynomial, BoundingBox(c,this->xmax(), this->ymin(), this->ymax(), this->zmin(), this->zmax()));
    /*  Update neighbors  */
    this->connect0(left, right);
    left->join0(right);
  } else if (v==1) {
    left = new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(), this->ymin(), c, this->zmin(), this->zmax()));
    right= new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(), c, this->ymax(), this->zmin(), this->zmax()));
    /*  Update neighbors  */
    this->connect1(left, right);
    left->join1(right);
  } else {//v==2
    left = new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(),this->ymin(),this->ymax(), this->zmin(), c));
    right= new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(),this->ymin(),this->ymax(), c, this->zmax()));
    /*  Update neighbors  */
    this->connect2(left, right);
    left->join2(right);
  } 
  //bcell3d_split(left,right,this,v,c);
  //  cout<<"Subdivide end "<<*left<<" "<<*right<<endl;
  tensor::split(left->m_polynomial, right->m_polynomial, v);
  //this->disconnect();
}

//--------------------------------------------------------------------
TMPL void 
SELF::subdivide(CellBase*& Left, CellBase*& Right, int v, double c) {
  SELF * left=0, * right=0;
  this->split(left,right,v,c);
  Left=left; Right=right;
}

//--------------------------------------------------------------------
TMPL bool 
SELF::insert_regular(Topology* t) {
  
  //  Topology3d* tpl3d= dynamic_cast<Topology3d*>(t);
  //  marching_cube::polygonise(*this, *this);
    
  //  if(this->m_idx == 0) 
  
  // marching_cube::polygonise(*t, *this); //NOT before topology_regular
  //  tpl3d->insert(F);
  
  
  //foreach(Point* p, F->points()) tpl3d->insert(p);
  return true;
}
//--------------------------------------------------------------------
TMPL void
SELF::polygonise(Topology3d* t) {
  marching_cube::polygonise(*t, *this); 
}
//--------------------------------------------------------------------
TMPL bool 
SELF::insert_singular(Topology* ) {
  return true;
}
//--------------------------------------------------------------------
TMPL bool 
SELF::topology_regular(Topology* t) {
  
  Topology3d* tpl3d= dynamic_cast<Topology3d*>(t);

   if(this->w_neighbors.size()>1)
     resolve_conflict(this->w_neighbors, this->xmin(),0 );

   if(this->e_neighbors.size()>1)
     resolve_conflict(this->e_neighbors, this->xmax(),0 );

    // foreach(Cell3d* cl, this->n_neighbors) 
    //   if(dynamic_cast<SELF*>(cl)->m_faces.size()>0)
    //     insert_bbx(tpl3d,cl,1,0);      
  
  if(this->n_neighbors.size()>1)
    resolve_conflict(this->n_neighbors, this->ymax(),1 );

   if(this->s_neighbors.size()>1)
     resolve_conflict(this->s_neighbors, this->ymin(),1 );

   if(this->f_neighbors.size()>1)
     resolve_conflict(this->f_neighbors, this->zmax(),2 );

   if(this->b_neighbors.size()>1)
     resolve_conflict(this->b_neighbors, this->zmin(),2 );

    //output topology
  foreach(Point* p, this->m_points)  tpl3d->insert(p);
  foreach(Face* f,  this->m_faces)   
  {
    tpl3d->insert(f);
    //if (f->size()>3) insert_bbx(tpl3d, (BoundingBox*)this);
  }

  return true;
}

//--------------------------------------------------------------------
TMPL void
SELF::resolve_conflict(Seq<Cell3d*> &neighb, double side, int cc)
{

  typedef Cell3d Cell3dRef;
    Seq<Point*> tmp;
    Seq<Edge*> edges;
    bool ccw;
    
    //if (0)
    foreach(Face* f1, this->m_faces) 
    { // discover wrong egde
      ccw= f1->is_ccw();
      tmp.clear();
      foreach(Point *p, *f1)
        if ( (*p)[cc]==side )  //( p->y()==this->ymax() )
          tmp<<p;
      if ( tmp.size()==2 )
      {
        Edge * tt=new Edge(tmp[0],tmp[1]);
        edges.clear();
        foreach(Cell3dRef* cl, neighb) //traverse neighbors 
          if(dynamic_cast<SELF*>(cl)->m_faces.size()>0)
          {
            foreach(Face* f0,dynamic_cast<SELF*>(cl)->m_faces)
            {// discover right edges
              tmp.clear();
              foreach(Point *p, *f0)
                //if ( p->y() == cl->ymin()  )
                if ( (*p)[cc] == side )
                  tmp<<p;
              if ( tmp.size()==2 )
              { //gather edge
                edges<< new Edge( tmp[0],tmp[1]);
              }
            }
          }
        if ( edges.size()>0 )
        {
          std::cout<<"Found big edge on plane "<<cc<<", "<<side<<std::endl;
          std::cout<<"adding "<<edges.size() << " edges." <<std::endl;
        }
        foreach( Edge* e, edges) // Split f1 by adding new edges
        {
          Face * ff=
            new Face( 
              e->destination(), 
              e->source(), 
              f1->not_on(tt) );
          if ( ff->is_ccw()!=ccw )
            ff->reverse();
          this->m_faces <<  ff;
        }
        m_faces.erase( m_faces.search(f1) );//delete old face
      }
    }
}

//--------------------------------------------------------------------
TMPL double
SELF::vertex_eval(unsigned sx, unsigned sy, unsigned sz) const {
  
  int s=0;
  s+= sx*(m_polynomial.rep().env.sz(0)-1)*m_polynomial.rep().env.st(0);
  s+= sy*(m_polynomial.rep().env.sz(1)-1)*m_polynomial.rep().env.st(1);
  s+= sz*(m_polynomial.rep().env.sz(2)-1)*m_polynomial.rep().env.st(2);
  return m_polynomial[s];
}

//--------------------------------------------------------------------
TMPL int
SELF::mc_index() const {

  double isolevel = 0;
  int CubeIndex=0;
  
  if (vertex_eval(0,1,0) <= isolevel) CubeIndex |= 1;
  if (vertex_eval(1,1,0) <= isolevel) CubeIndex |= 2;
  if (vertex_eval(1,0,0) <= isolevel) CubeIndex |= 4;
  if (vertex_eval(0,0,0) <= isolevel) CubeIndex |= 8;
  if (vertex_eval(0,1,1) <= isolevel) CubeIndex |= 16;
  if (vertex_eval(1,1,1) <= isolevel) CubeIndex |= 32;
  if (vertex_eval(1,0,1) <= isolevel) CubeIndex |= 64;
  if (vertex_eval(0,0,1) <= isolevel) CubeIndex |= 128;
  
  return CubeIndex;
}
//--------------------------------------------------------------------
TMPL bool
SELF::edge_point(Point** CELL, int cube_index) const {
  
  BoundingBox* bx = (BoundingBox*)this;
  Seq<Point*>  F;
  Polynomial bk, ft;

  tensor::face(bk, equation(), 2, 0);
  if (cube_index & 1) {
    Solver::edge_point(F, bk, Solver::north_back_edge, *bx);    
    if(F.size()==0) { 
      std::cout<<"Problem in MC0"<<std::endl;
      return false;
    } else
      CELL[0] = F[0]; 
    F.resize(0);
  }
  if (cube_index & 2) {
    Solver::edge_point(F, bk, Solver::east_back_edge , *bx);
    if(F.size()==0) {  
      std::cout<<"Problem in MC1"<<std::endl;
      return false;
    } else
      CELL[1] = F[0];  
    F.resize(0);
  }
  if (cube_index & 4) {
    Solver::edge_point(F, bk, Solver::south_back_edge, *bx);
    if(F.size()==0) {  
      std::cout<<"Problem in MC2"<<std::endl;
      return false;
    } else
      CELL[2] = F[0]; 
    F.resize(0);
  }
  if (cube_index & 8) {
    Solver::edge_point(F, bk, Solver::west_back_edge , *bx);
    if(F.size()==0) { 
      std::cout<<"Problem in MC3"<<std::endl;
    } else
      CELL[3] = F[0]; 
    F.resize(0);
  }

  tensor::face(ft,equation(), 2, 1);

  if (cube_index & 16) {
    Solver::edge_point(F, ft, Solver::north_front_edge, *bx);
    if(F.size()==0) { 
      std::cout<<"Problem in MC4 "<<std::endl;
      return false;
      } else
      CELL[4] = F[0]; 
    F.resize(0);
  }
  
  if (cube_index & 32) {
    Solver::edge_point(F, ft, Solver::east_front_edge , *bx);
    if(F.size()==0) { 
      std::cout<<"Problem in MC5"<<std::endl;
      return false;
    } else
      CELL[5] = F[0]; 
    F.resize(0);
  }
  
  if (cube_index & 64) {
    Solver::edge_point(F, ft, Solver::south_front_edge, *bx);
    if(F.size()==0) {
      std::cout<<"Problem in MC6"<<std::endl;
      return false;
    } else
      CELL[6] = F[0]; 
    F.resize(0);
  }
  
  if (cube_index & 128) {
    Solver::edge_point(F, ft, Solver::west_front_edge , *bx);
    if(F.size()==0)  {
      std::cout<<"Problem in MC7"<<std::endl;
      return false;
    } else
      CELL[7] = F[0]; 
    F.resize(0);
  }
  
  tensor::face(ft, equation(), 1, 1);
  if (cube_index & 256) {
    Solver::edge_point(F, ft, Solver::north_west_edge, *bx);
    if(F.size()==0)  {
      std::cout<<"Problem in MC8"<<std::endl;
      return false;
    } else
      CELL[8] = F[0]; 
    F.resize(0);
  } 
  
  if (cube_index & 512) {
    Solver::edge_point(F, ft, Solver::north_east_edge, *bx);
    if(F.size()==0)  {
      std::cout<<"Problem in MC9"<<std::endl;
      return false;
    } else
      CELL[9] = F[0];
    F.resize(0);
  } 
  
  tensor::face(bk, equation(), 1, 0);
  
  if (cube_index & 1024) {
    Solver::edge_point(F, bk, Solver::south_east_edge, *bx);
    if(F.size()==0)  {
      std::cout<<"Problem in MC10"<<std::endl;
      return false;
    } else
      CELL[10] = F[0]; 
    F.resize(0);
    } 
  
  if (cube_index & 2048) {
    Solver::edge_point(F, bk, Solver::south_west_edge, *bx);
    if(F.size()==0)  {
      std::cout<<"Problem in MC11"<<std::endl;
      return false;
    } else
      CELL[11] = F[0]; 
    F.resize(0);
  } 

  return true;
}

//====================================================================
} ; // namespace shape
} ; // namespace mmx
# undef TMPL
# undef TMPL1
# undef Point
# undef Edge
# undef AlgebraicSurface
# undef SELF 
# endif // SHAPE_IMPLICITCELL_H

---