synaps/topology/Cell3dN.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_cell3dN_H
#define synaps_topology_cell3dN_H

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

#include <synaps/mpol/MPol.h>
#include <synaps/mpol/MPolDse.h>
#include <synaps/topology/box3d.h>
#include <synaps/topology/inter_points.h>
#include <synaps/topology/point_graph.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
//--------------------------------------------------------------------

//====================================================================
// Cell3dN: the structure of the cell that is tailored for building topological graph of
//               the curve defined by two intersection implicit surfaces
template < int N, 
           class C = double,
           class MP = MPolDse<C, mpoldse::bernstein<C> > >
struct Cell3dN
{
  typedef C                          coeff_t;
  typedef topology::box3d<coeff_t>   box_t;
  typedef MP                         MPOL_T;
  //  typedef topology::point<C>         point_t;
  typedef Cell3dN<N,C,MP>  cell_t;
  typedef Cell3dN<N,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;
  MPOL_T      m_F[N]; // surface functions
  MPOL_T*     m_T;//=M_F+2; // tangent fields
  
  bool        m_aFaceComputed[6]; // flag indicating wether the intersection has been computed for a facet
  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 //
  
  Cell3dN();
  Cell3dN(const coeff_t& a0, const coeff_t& b0, const coeff_t& a1, const coeff_t& b1, const coeff_t& a2, const coeff_t& b2);
  Cell3dN(const cell_t& c); // copy constructor
  
  // operator overrides //
  
  // assign operator
  self_t&  operator = (const cell_t& c);
  
  // helpers  //

  int nbpol() const {return N;}
  // 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 <int N, class C, class MP >
Cell3dN<N, C, MP>::Cell3dN() :
  m_T(m_F+2), m_nDepth(0), m_nRegularIdx(-2)
{
    setBox(0,0,0,0,0,0);
    memset(m_aFaceComputed, false, sizeof(bool)*6);
}


template <int N, class C, class MP >
Cell3dN<N, C, MP>::Cell3dN(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_T(m_F+2), 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_aFaceComputed, false, sizeof(bool)*6);
}

template <int N, class C, class MP >
Cell3dN<N, C, MP>::Cell3dN(const cell_t& c) :
  m_T(m_F+2), m_nDepth(c.m_nDepth), m_nRegularIdx(c.m_nRegularIdx)
{
    assign(c);
}


// operator overrides //

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

// helpers //
//--------------------------------------------------------------------
template < int N, class C, class MP >
void Cell3dN<N, C, MP>::assign(const cell_t& c)
{
    m_cBox = c.m_cBox;
    for (int i = 0; i < N; i++)
        m_F[i] = c.m_F[i];
    m_T = m_F+2;
//     for (int i = 0; i < 3; i++)
//         m_T[i] = c.m_T[i];

    m_nDepth = c.m_nDepth;
    m_nRegularIdx = c.m_nRegularIdx;
    memcpy(m_aFaceComputed, c.m_aFaceComputed, sizeof(bool)*6);
    m_aInterList = c.m_aInterList;
}
//--------------------------------------------------------------------
template < int N, class C, class MP >
void Cell3dN<N, 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 < int N, class C, class MP >
void Cell3dN<N, 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 < int N, class C, class MP>
void Cell3dN<N, 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 <int N, class C, class MP>
bool Cell3dN<N,C,MP>::containsVariety(unsigned i) const
{
  assert(i<N);
  bool hasroot;

  const C* coeff = this->m_F[i].begin();
  C lastC = coeff[0];
  hasroot = FALSE;
  
  for (int j = 1; j < this->m_F[i].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;
        }
    }
  
  if (!hasroot) // if either surface (0-levelset) does not cross this cell, returns FALSE
    return FALSE;
  return TRUE;
}
//--------------------------------------------------------------------
// test whether its zero-levelset crosses the current cell or not
template <int N, class C, class MP>
template <class POINT>
bool Cell3dN<N,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 <int N, class C, class MP >
std::ostream & operator<<( std::ostream & os, const Cell3dN<N, C, MP>& c )
{
    os  << "Coord: " << c.m_cBox << std::endl
        << "F: " << c.m_F[0] << std::endl
        << "G: " << c.m_F[1] << std::endl
        << "t1: " << c.m_T[0] << std::endl
        << "t2: " << c.m_T[1] << std::endl
        << "t3: " << c.m_T[2] << std::endl;
    return os;
}
//====================================================================
__END_NAMESPACE_SYNAPS
/********************************************************************/
#endif //synaps_topology_cell3dN_H