synaps/topology/Cell3d.h

/***************************************************************************
 *   Copyright (C) 2005 by Chen Liang                                      *
 *   cliang@cs.hku.hk                                                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifndef synaps_topology_cell3d_H
#define synaps_topology_cell3d_H

#include <vector>
#include <iostream>
#include <assert.h>
#include <vector>

#include <synaps/mpol/MPol.h>
#include <synaps/mpol/MPolDse.h>
#include <synaps/topology/box3d.h>
#include <synaps/topology/inter_points.h>


// custom macros
#define VARIDX(poly, idx) (idx)
//#define VARIDX(poly, idx) (poly.rep().env.vr(idx))


#define TRUE true
#define FALSE false

#ifdef SYNAPS_SBD3D_USEGMP
#ifdef SYNAPS_WITH_GMP
#define _is_singular(x) ((x) == QQ("0"))
#define _is_verysmall(x) (abs(x)<QQ("0.001"))
#else
#define _is_singular(x) ((x) < 1e-6 && (x) > -1e-6)
#define _is_verysmall(x) (std::abs(x)<1e-3)
#endif //end: SYNAPS_WITH_GMP
#else
#define _is_singular(x) ((x) < 1e-6 && (x) > -1e-6)
#define _is_verysmall(x) (std::abs(x)<1e-3)
#endif //end: SYNAPS_SBD3D_USEGMP

#ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
// random factor
#include <stdlib.h> // for rnd
#endif


//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//====================================================================
// Cell3d: the structure of the cell that is tailored for building topological graph of
//               the curve defined by two intersection implicit surfaces
template < class C = double,
           class MP = MPolDse<C, mpoldse::bernstein<C> > >
struct Cell3d
{
  typedef C                          coeff_t;
  typedef topology::box3d<coeff_t>   box_t;
  typedef MP                         MPOL_T;
  //  typedef topology::point<C>         point_t;
  typedef Cell3d<C,MP>  cell_t;
  typedef Cell3d<C,MP>  self_t;
  
  typedef intersection3d<C>        intersection_t;
  typedef Seq<intersection3d<C>*>  intersectionlist_t;
  typedef shared_object<intersection_t*>  intersection_ptr;
  typedef Seq<intersection_ptr>   intersection_ptr_list;
  typedef topology::point_graph<C, intersection_ptr_list> pointgraph_int_t;
  
  // data members //

  box_t       m_cBox;
  std::vector<MPOL_T>  m_F; // surface functions
  unsigned             m_nT; // index of tangent fields
  
  bool                m_FaceComputed[6]; // flag indicating wether the intersection has been computed for a facet
  intersectionlist_t  m_FaceLstPtIdx[6];  // list of points on the corresponding face, including the edges.

  intersectionlist_t  m_EdgeLstPtIdx[12];  // list of points on the corresponding face, including the edges.


  int                 m_nDepth;      // the depth of this cell in the octree
  int                 m_nRegularIdx; // -2 = not computed, -1 = singular, 0-k = regular w.r.t.
  intersectionlist_t  m_aInterList; // a list of intersections already computed for this cell
  
  // constructors //
  
  //  Cell3d();
  Cell3d(unsigned N);
  Cell3d(unsigned N, const coeff_t& a0, const coeff_t& b0, const coeff_t& a1, const coeff_t& b1, const coeff_t& a2, const coeff_t& b2);
  Cell3d(const MPOL_T& p, const coeff_t& a0, const coeff_t& b0, const coeff_t& a1, const coeff_t& b1, const coeff_t& a2, const coeff_t& b2);
  //  Cell3d(const coeff_t& a0, const coeff_t& b0, const coeff_t& a1, const coeff_t& b1, const coeff_t& a2, const coeff_t& b2);  
  Cell3d(const cell_t& c); // copy constructor
  
  // operator overrides //
  
  // assign operator
  self_t&  operator = (const cell_t& c);
  
  // helpers  //

  int nbpol() const {return m_F.size();}
  void push_back(const MPOL_T& p) { m_F.push_back(p); }
  // this function copies the content of cell c to this cell, used by operator = and copy constructor
  void     assign(const cell_t&);

  // merge two lists of type intersectionlist_t
  void     mergelist(intersectionlist_t& dstlist, const intersectionlist_t& srclist);

  // merge two lists of type intersectionlist_t with points inside the cell
  void     merge_inside(intersectionlist_t& dstlist, const intersectionlist_t& srclist);
  
  // member functions //
  
  // this function set the bounding box (or coordinate represented by this cell)
  virtual void setBox(const C& a0, const C& b0, 
                      const C& a1, const C& b1, 
                      const C& a2, const C& b2);

  inline virtual bool testsize(coeff_t epsilon)
  {
    if (m_cBox.x1-m_cBox.x0 > epsilon || m_cBox.z1-m_cBox.z0 > epsilon || m_cBox.z1-m_cBox.z0 > epsilon)
      {
        return FALSE;
      }
    return TRUE;
  }

  bool containsVariety(unsigned i) const;
  template<class POINT>
  bool contains_point(const POINT& pt) const;
};
//====================================================================
// constructors //

// template <class C, class MP >
// Cell3d<C, MP>::Cell3d() :
//   m_nT(1), m_nDepth(0), m_nRegularIdx(-2)
// {
//     setBox(0,0,0,0,0,0);
//     memset(m_FaceComputed, false, sizeof(bool)*6);
// }


template <class C, class MP >
Cell3d<C, MP>::Cell3d(unsigned N) :
  m_F(N), m_nT(1), m_nDepth(0), m_nRegularIdx(-2)
{
    setBox(0,0,0,0,0,0);
    memset(m_FaceComputed, false, sizeof(bool)*6);
}


// template <class C, class MP >
// Cell3d<C, MP>::Cell3d(const coeff_t& a0, const coeff_t& b0, 
//                    const coeff_t& a1, const coeff_t& b1, 
//                    const coeff_t& a2, const coeff_t& b2) :
//   m_nT(1), m_nDepth(0), m_nRegularIdx(-2)
// {
//     // m_cBox = [a0,b0]X[a1,b1]X[a2,b2]
//     setBox(a0,b0,a1,b1,a2,b2);
//     memset(m_FaceComputed, false, sizeof(bool)*6);
// }

template <class C, class MP >
Cell3d<C, MP>::Cell3d(const MP & p,
                      const coeff_t& a0, const coeff_t& b0, 
                      const coeff_t& a1, const coeff_t& b1, 
                      const coeff_t& a2, const coeff_t& b2) :
  m_nT(1), m_nDepth(0), m_nRegularIdx(-2)
{
    // m_cBox = [a0,b0]X[a1,b1]X[a2,b2]
    setBox(a0,b0,a1,b1,a2,b2);
    memset(m_FaceComputed, false, sizeof(bool)*6);
    m_F.push_back(p);
}

template <class C, class MP >
Cell3d<C, MP>::Cell3d(unsigned N,
                      const coeff_t& a0, const coeff_t& b0, 
                      const coeff_t& a1, const coeff_t& b1, 
                      const coeff_t& a2, const coeff_t& b2) :
  m_F(N), m_nT(1), m_nDepth(0), m_nRegularIdx(-2)
{
    // m_cBox = [a0,b0]X[a1,b1]X[a2,b2]
    setBox(a0,b0,a1,b1,a2,b2);
    memset(m_FaceComputed, false, sizeof(bool)*6);
}

template <class C, class MP >
Cell3d<C, MP>::Cell3d(const cell_t& c) :
  m_F(c.nbpol()), m_nT(1), m_nDepth(c.m_nDepth), m_nRegularIdx(c.m_nRegularIdx)
{
    assign(c);
}


// operator overrides //

template <class C, class MP >
inline Cell3d<C, MP>& Cell3d<C, MP>::operator = (const cell_t& c)
{
    assign(c); return *this;
}

// helpers //
//--------------------------------------------------------------------
template < class C, class MP >
void Cell3d<C, MP>::assign(const cell_t& c)
{
    m_cBox = c.m_cBox;
    m_F = c.m_F;
    m_nT = c.m_nT;
    m_nDepth = c.m_nDepth;
    m_nRegularIdx = c.m_nRegularIdx;
    memcpy(m_FaceComputed, c.m_FaceComputed, sizeof(bool)*6);
    m_aInterList = c.m_aInterList;
}
//--------------------------------------------------------------------
template <class C, class MP >
void Cell3d<C, MP>::mergelist(intersectionlist_t& dstlist, 
                              const intersectionlist_t& srclist)
{
    for (int v = 0; v < (int)srclist.size(); v++)
        dstlist.push_back(srclist[v]);
}

//--------------------------------------------------------------------
template <class C, class MP >
void Cell3d<C, MP>::merge_inside(intersectionlist_t& dstlist, 
                                 const intersectionlist_t& srclist)
{
  PRINT_DEBUG("box "<<m_cBox);
    for (int v = 0; v < (int)srclist.size(); v++)
      if(contains_point(srclist[v]->pnt))
        {
          PRINT_DEBUG(" "<<srclist[v]->pnt);
          dstlist.push_back(srclist[v]);
        }
}
//--------------------------------------------------------------------
template <class C, class MP>
void Cell3d<C, MP>::setBox(const coeff_t& a0, const coeff_t& b0, const coeff_t& a1, const coeff_t& b1, const coeff_t& a2, const coeff_t& b2)
{
    m_cBox = box_t(a0,b0,a1,b1,a2,b2);
}
//--------------------------------------------------------------------
// test whether its zero-levelset crosses the current cell or not
template <class C, class MP>
bool Cell3d<C,MP>::containsVariety(unsigned i) const
{
  assert(i<nbpol());
  return has_sign_change(this->m_F[i]);
}
//--------------------------------------------------------------------
// test whether its zero-levelset crosses the current cell or not
template <class C, class MP>
template <class POINT>
bool Cell3d<C,MP>::contains_point(const POINT& pt) const
{
  return ( pt[0] >= m_cBox.x0 && pt[0] <= m_cBox.x1 && 
           pt[1] >= m_cBox.y0 && pt[1] <= m_cBox.y1 && 
           pt[2] >= m_cBox.z0 && pt[2] <= m_cBox.z1);
}

// non-member functions //

template <class C, class MP >
std::ostream & operator<<( std::ostream & os, const Cell3d<C, MP>& c )
{
  os  << "Coord: " << c.m_cBox << std::endl;
  for(unsigned i=0;i<c.nbpol();i++)
    os << " " <<c.m_F[i] << std::endl;
  
  return os;
}
//--------------------------------------------------------------------
template <class MP>
bool has_sign_change(const MP& pol)
{
  typedef typename MP::coeff_t C;
  bool hasroot;
  
  const C* coeff = pol.begin();
  C lastC = coeff[0];
  hasroot = FALSE;
  
  for (int j = 1; j < pol.rep().env.sz(); j++)
    {
      C cirC = coeff[j];
      if (_is_singular(cirC)) // treat it as 0
        {
          continue;
        }
      else
        {
          if (lastC*cirC < 0) // a sign change of the coefficient is detected
            {
              hasroot = TRUE;
              break;
            }
          lastC = cirC;
        }
    }
  
  return hasroot;
}

template<class PT, class C, class MP>
bool is_on_face(const PT& pt, Cell3d<C,MP>* c, unsigned i, unsigned j)
{
  if(i==0)
    if (j==0)
      return( pt[0]== c->m_cBox.x0);
    else
      return( pt[0]== c->m_cBox.x1);
  else if(i==1)
    if(j==0)
      return( pt[1]== c->m_cBox.y0);
    else
      return( pt[1]== c->m_cBox.y1);
  else if(i==2) 
    if(j==0)
      return( pt[2]== c->m_cBox.z0);
    else  //if(i==2 && j==1)
      return( pt[2]== c->m_cBox.z1);
  return false;
}

template<class PT, class C, class MP>
bool is_on_border(const PT& pt, Cell3d<C,MP>* c, unsigned i)
{
  switch(i)
    {
    case 0:
      return((pt[1]== c->m_cBox.y0) || (pt[1]== c->m_cBox.y1) || (pt[2]== c->m_cBox.z0) || (pt[2]== c->m_cBox.z1));
    case 1:
      return((pt[0]== c->m_cBox.x0) || (pt[0]== c->m_cBox.x1) || (pt[2]== c->m_cBox.z0) || (pt[2]== c->m_cBox.z1));
    case 2:
      return((pt[0]== c->m_cBox.x0) || (pt[0]== c->m_cBox.x1) || (pt[1]== c->m_cBox.y0) || (pt[1]== c->m_cBox.y1));
    }
}
template<class PT, class C, class MP>
bool is_edge(const PT& pt, Cell3d<C,MP>* c, unsigned i)
{
  switch(i)
    {
    case 0:
      return((pt[1]== c->m_cBox.y0) || (pt[1]== c->m_cBox.y1) || (pt[2]== c->m_cBox.z0) || (pt[2]== c->m_cBox.z1));
    case 1:
      return((pt[0]== c->m_cBox.x0) || (pt[0]== c->m_cBox.x1) || (pt[2]== c->m_cBox.z0) || (pt[2]== c->m_cBox.z1));
    case 2:
      return((pt[0]== c->m_cBox.x0) || (pt[0]== c->m_cBox.x1) || (pt[1]== c->m_cBox.y0) || (pt[1]== c->m_cBox.y1));
    }
}



// template<class PT, class C, class MP>
// OCTREE::FACE_TYPE face_of(const PT& pt, const Cell3d<C,MP> & c)
// {
//   if (pt[0]== c.m_cBox.x0)
//     return B;
//   else if (pt[0]== c.m_cBox.x1)
//     return F;
//   else if (pt[1]== c.m_cBox.y0)
//     return W;
//   else if (pt[1]== c.m_cBox.y1)
//     return E;
//   else if (pt[2]== c.m_cBox.z0)
//     return S;
//   else 
//     return N:
// }
//====================================================================
__END_NAMESPACE_SYNAPS
/********************************************************************/
#endif //synaps_topology_cell3d_H