CGAL users discussion list

[Laurent Saboret <>]

Mariette Yvinec wrote:
> Qianqian Fang wrote:
> > hi
> >
> > In the past, I asked in this list about removing self-intersecting
> > elements from triangular surface. I was recommended to use
> > the surface mesher avoid repairing a bad mesh as it guarantees to
> > have no self-intersecting elements.
> >
> > Now, I need to process a surface generated from other software,
> > unfortunately, this mesh contains self-intersecting elements. So, I
> > still need a program to fix this surface (and resample it to my
> > desired density).
> >
> > Reading the documents of the Surface Mesher, I have the impression
> > that it can be applied for this situation. The only difference
> > is that the surface is not defined by gray-scale image, or an implicit
> > function, but by a triangular surface mesh.
> >
> > Indeed, I can even define the surface by an implicit function, which
> > returns the distance from a Point_3 to a Polyhedron_3, similar to
> > the sphere case in the manual. Does this sound like a viable
> > approach to you?
> >
> > If it may work, what I need is to find a subroutine to calculate the
> > distance. I only found a squared_distance() for Point_3
> > to Plane_3 in include/CGAL/squared_distance_3_2.h. I am wondering
> > if a similar function to calculate Point3 to a polyhedron exists?
> > if it does, can anyone show me an example how to use it?
> >
> > thank you
> >
> > Qianqian
>
> The Surface mesher will be able to handle such a situation in a
> a near future.
> What is missing actually is the Polyhedral_surface_traits that can handle
> surface defined as a triangle soup. Such a traits is going to be provided
> in the next public release.
>
> In the current release, you can find such a traits hidden in the polyhedron
> demo, as Pierre said.  However you will have to hack a bit the code, ti 
> bypass the construction of the Polyhedron_3 which does not work with a self 
> intersecting triangular surface.

Hi Qianqian Fang,

I attached to this post the AABB_polyhedral_oracle class mentioned by Mariette above, that 
can remesh a triangles soup. Triangles must be convertible to type CGAL::Triangle_3.

Best regards,
--
Laurent Saboret
INRIA Sophia-Antipolis


#ifndef AABB_POLYHEDRAL_ORACLE_H
#define AABB_POLYHEDRAL_ORACLE_H

#include <boost/static_warning.hpp>
#include <utility>
#include <CGAL/iterator.h>

#include "AABB_tree.h"


// Modified version of CGAL::AABB_polyhedral_oracle to be independent of a 
polyhedron class.
template <class TriangleKernel,           ///< Triangles' kernel
          class TriangleIterator,         ///< Value type must be Triangle_3
          class AABBTree_kernel = CGAL::Simple_cartesian<typename 
TriangleKernel::FT>
                                          ///< Local kernel for fast bounding 
box intersections
>                                          
//template <class Polyhedron, class TriangleKernel, class AABBTree_kernel>
class AABB_polyhedral_oracle 
//: public Polyhedron
{
    typedef AABB_polyhedral_oracle Self;
    
public:
    //typedef typename TriangleKernel::FT FT;
    //typedef typename TriangleKernel::Ray_3 Ray_3;
    //typedef typename TriangleKernel::Line_3 Line_3;
    typedef typename TriangleKernel::Point_3 Point_3;
    //typedef typename TriangleKernel::Segment_3 Segment_3;

    typedef Self Surface_mesher_traits_3;
    typedef Point_3 Intersection_point;
    typedef Self Surface_3;

    // AABB tree
    typedef AABB_tree<TriangleKernel, TriangleIterator, AABBTree_kernel> Tree;
    typedef typename Tree::Point_with_triangle Point_with_triangle;

private:
    Tree *m_pTree;

public:
    Tree* tree() const{ return m_pTree; }

public:
    // constructor
    AABB_polyhedral_oracle()
    {
        m_pTree = NULL;
    }
    AABB_polyhedral_oracle(Tree *pTree)
    {
        m_pTree = pTree;
    }
    //AABB_polyhedral_oracle(const AABB_polyhedral_oracle& oracle)
    //{
    //    m_pTree = oracle.tree();
    //}

    //class Intersect_3;
    //friend class Intersect_3;

    class Intersect_3 {
        const Self& self;
    public:
        Intersect_3(const Self& self) : self(self)
        {
        }

        template<class SegmentType>    ///< tested with Segment_3 and Ray_3
        CGAL::Object operator()(const Surface_3& surface, const SegmentType& 
segment) const
        {
            Point_with_triangle pwh;
            if(surface.tree()->first_intersection(segment,pwh))
                return CGAL::make_object(pwh.first);
            else
                return CGAL::Object();
        }

        //Object operator()(const Surface_3& surface, const Ray_3& ray) const
        //{
        //    Converter convert;
        //    BConverter bconvert;
        //    Point_with_triangle pwh;
        //    if(surface.tree()->first_intersection(convert(ray),pwh))
        //        return make_object(bconvert(pwh.first));
        //    else
        //        return Object();
        //}

        //Object operator()(const Surface_3& surface, const Line_3& line) 
const
        //{
        //    Converter convert;
        //    BConverter bconvert;
        //    Point_with_triangle pwh;
        //    if(surface.tree()->first_intersection(convert(line),pwh))
        //        return make_object(bconvert(pwh.first));
        //    else
        //        return Object();
        //}
        
    }; // class Intersect_3

    Intersect_3 intersect_3_object() const
    {
        return Intersect_3(*this);
    }

    //class Construct_initial_points;
    //friend class Construct_initial_points;

    class Construct_initial_points
    {
        const Self& self;
    public:
        Construct_initial_points(const Self& self) : self(self)
        {
        }

        template <typename OutputIteratorPoints>
        OutputIteratorPoints operator() (const Surface_3& surface, 
                                         OutputIteratorPoints out, 
                                         int n) const 
        {
            // TODO (with visitor)
            // std::cout << "construct initial point set" << std::endl;
            // *out++= p;
            return out;
        }
        
    }; // class Construct_initial_points

    Construct_initial_points construct_initial_points_object() const
    {
        return Construct_initial_points(*this);
    }

    template <class P>
    bool is_in_volume(const Surface_3& surface, const P& p)
    {
        std::cout << "DEBUG: call is in volume" << std::endl;
        return true;
    }
    
}; // end class AABB_polyhedral_oracle


#endif // AABB_POLYHEDRAL_ORACLE_H
// Author(s)     :  Camille Wormser, Jane Tournois, Pierre Alliez, Laurent 
Saboret

#ifndef AABB_TREE_H_INCLUDED
#define AABB_TREE_H_INCLUDED

#include "AABB_node.h"

#include <CGAL/Simple_cartesian.h>

#include <stack>


// Modified version of CGAL::AABB_tree to be independent of a polyhedron 
class.
template <class TriangleKernel,           ///< Triangles' kernel
          class TriangleIterator,         ///< Value type must be Triangle_3
          class Kernel = CGAL::Simple_cartesian<typename TriangleKernel::FT>
                                          ///< Local kernel for fast bounding 
box intersections
>                                          
class AABB_tree
{
private:

  // Nodes of the tree
  typedef AABB_node<TriangleKernel, 
                    TriangleIterator,  
                    typename std::vector<TriangleIterator>::iterator,
                    Kernel> Node;

public:

  // Local kernel object types
  typedef typename Kernel::FT FT;
  typedef typename Kernel::Point_3 Point;

  // Global kernel object types
  typedef typename TriangleKernel::Triangle_3 Triangle;
  typedef typename TriangleKernel::Point_3 Triangle_point;
  typedef std::pair<Triangle_point, TriangleIterator> Point_with_triangle;

private:

  // set of triangle pointers (global to the tree)
  std::vector<TriangleIterator> m_data;

  // root node
  Node* m_root;

public:

  // life cycle
  AABB_tree() : m_root(NULL) {}
  AABB_tree(TriangleIterator first, TriangleIterator beyond)
  {
    m_root = NULL;
    build_faces(first, beyond);
  }

  ~AABB_tree()
  {
    cleanup();
  }

  void cleanup()
  {
    m_data.clear();
    delete m_root; m_root = NULL;
  }

  // build tree
  bool build_faces(TriangleIterator first, TriangleIterator beyond)
  {
    cleanup();
    set_face_data(first, beyond);
    if(!empty())
    {
      m_root = new Node(m_data.begin(), m_data.end());
      return true;
    }
    return false;
  }

private:

  void set_face_data(TriangleIterator first, TriangleIterator beyond)
  {
    unsigned int nbf = std::distance(first, beyond);
    m_data.reserve(nbf);
    for (TriangleIterator triangle = first; triangle != beyond; triangle++)
      m_data.push_back(triangle);
  }

public:

  bool empty()
  {
    return m_data.size() < 2; // TODO: change this requirement to < 1
  }

  // --------------------RAY/SEGMENT ORACLES----------------------//

  template<class SegmentType>    ///< tested with Segment_3 and Ray_3
  bool first_intersection(const SegmentType& seg, Point_with_triangle& pwh) 
const
  {
    std::stack<Node*> nodes;
    nodes.push(m_root);
    while(!nodes.empty())
    {
      Node *pNode = nodes.top();
      nodes.pop();

      // Test against primitive if leaf, against box if not
      if(pNode->do_intersect(seg)) 
      {
        if(pNode->is_leaf())
        {
          // found intersecting triangle
          Triangle t(*pNode->primitive());
          CGAL::Object inter = CGAL::intersection(t.supporting_plane(),seg);
          Triangle_point result;
          bool ok = CGAL::assign(result, inter);
          assert(ok);
          pwh = Point_with_triangle(result, pNode->primitive());
          return true; 
        }
        else // recurse (test against bbox only)
        {
          nodes.push(pNode->left_child());
          nodes.push(pNode->right_child());
        }
      }
    }
    return false;
  }


  template<class SegmentType,    ///< tested with Segment_3 and Ray_3
           class OutputIterator> ///< value type must be TriangleIterator
  OutputIterator find_intersections(const SegmentType& seg, OutputIterator 
output) const
  {
    std::stack<Node*> nodes;
    nodes.push(m_root);
    while(!nodes.empty())
    {
      Node *pNode = nodes.top();
      nodes.pop();

      // Test against primitive if leaf, against box if not
      if(pNode->do_intersect(seg)) 
      {
        if(pNode->is_leaf())
        {
          *output = pNode->primitive();
          output++; 
        }
        else // recurse (test against bbox only)
        {
          nodes.push(pNode->left_child());
          nodes.push(pNode->right_child());
        }
      }
    }
    return output;
  }

  template<class SegmentType>     ///< tested with Segment_3 and Ray_3
  unsigned int nb_intersections(const SegmentType& seg) const
  {
    std::vector<TriangleIterator> intersections;
    find_intersections(seg, std::back_inserter(intersections));
    return intersections.size();
  }

}; // class AABB_tree


#endif // AABB_TREE_H_INCLUDED
// Author(s)     :  Camille Wormser, Jane Tournois, Pierre Alliez, Laurent 
Saboret

#ifndef AABB_NODE_H_INCLUDED
#define AABB_NODE_H_INCLUDED

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/intersections.h>

#include <limits>
#include <vector>
#include <cassert>


// Modified version of CGAL::AABB_node to be independent of a polyhedron 
class.
template <class TriangleKernel,           ///< Triangles' kernel
          class TriangleIterator,         ///< Value type must be Triangle_3
          class TriangleIteratorIterator, ///< Value type must be 
TriangleIterator
          class Kernel = CGAL::Simple_cartesian<typename TriangleKernel::FT>
                                          ///< Local kernel for fast bounding 
box intersections
>                                          
class AABB_node
{
public:

  // Bounding box with (double) floating point arithmetic
  typedef CGAL::Iso_cuboid_3<Kernel> Bbox;
  
  // Local kernel object types
  typedef typename Kernel::FT FT;
  typedef typename Kernel::Point_3 Point;

  // Global kernel object types
  typedef typename TriangleKernel::Triangle_3 Triangle;
  typedef typename TriangleKernel::Point_3 Triangle_point;
  
private:

  // Converter of number types, points and vectors
  typedef CGAL::Cartesian_converter<TriangleKernel, Kernel> Converter; 
  
private:

  // bounding box
  Bbox m_bbox;

  // if leaf: data
  TriangleIterator m_primitive;

  // else: children data as iterator range (stored only in the tree)...
  TriangleIteratorIterator m_begin, m_end;
  // ...and children nodes.
  AABB_node *m_left_child;
  AABB_node *m_right_child;

public:

  // life cycle
  AABB_node(TriangleIteratorIterator a, TriangleIteratorIterator b)
  {
    m_left_child = m_right_child = NULL;
    expand(a, b);
  }

  ~AABB_node()
  {
    cleanup_recurse();
  }

  // access to primitive
  TriangleIterator primitive()  { return m_primitive; }

  // access to children
  AABB_node* left_child()  { return m_left_child; }
  AABB_node* right_child() { return m_right_child; }

  bool is_leaf()
  {
    return m_left_child == NULL && m_right_child == NULL;
  }

  // deletes the whole subtree rooted at this node (except this node, of 
course).
  void cleanup_recurse()
  {
    delete m_left_child;
    delete m_right_child;
    m_left_child = m_right_child = NULL;
  }

private:

  // compute bbox for iterator range of input primitives
  void compute_bbox()
  {
    assert(m_begin != m_end);
    
    FT xmin,xmax,ymin,ymax,zmin,zmax;
    std::vector<TriangleIterator>::iterator it = m_begin;

    // init with any point in container
    TriangleIterator triangle = *it;
    const Triangle_point& any_point = (*triangle)[0];
    xmin = xmax = any_point.x();
    ymin = ymax = any_point.y();
    zmin = zmax = any_point.z();

    // compute bbox
    for(; it != m_end; it++)
    {
      TriangleIterator triangle = *it;
      xyz_min(triangle,xmin,ymin,zmin);
      xyz_max(triangle,xmax,ymax,zmax);
    }
    
    // Enlarge it a little to avoid missing intersections
    // due to floating point computation errors.
    FT epsilon = FT(100)*std::numeric_limits<FT>::epsilon();
    Point p(xmin-epsilon, ymin-epsilon, zmin-epsilon);
    Point q(xmax+epsilon, ymax+epsilon, zmax+epsilon);
    m_bbox = Bbox(p,q);
  }

  void xyz_min(TriangleIterator triangle,
               FT& xmin, FT& ymin, FT& zmin)
  {
    Triangle_point a = (*triangle)[0];
    Triangle_point b = (*triangle)[1];
    Triangle_point c = (*triangle)[2];
    xmin = (std::min)(xmin,(std::min)(a.x(),(std::min)(b.x(),c.x())));
    ymin = (std::min)(ymin,(std::min)(a.y(),(std::min)(b.y(),c.y())));
    zmin = (std::min)(zmin,(std::min)(a.z(),(std::min)(b.z(),c.z())));
  }

  void xyz_max(TriangleIterator triangle,
               FT& xmax, FT& ymax, FT& zmax)
  {
    Triangle_point a = (*triangle)[0];
    Triangle_point b = (*triangle)[1];
    Triangle_point c = (*triangle)[2];
    xmax = (std::max)(xmax,(std::max)(a.x(),(std::max)(b.x(),c.x())));
    ymax = (std::max)(ymax,(std::max)(a.y(),(std::max)(b.y(),c.y())));
    zmax = (std::max)(zmax,(std::max)(a.z(),(std::max)(b.z(),c.z())));
  }

  // builds the tree by recursive expansion.
  // [a,b[ is the range of primitives to be added to the tree.
  void expand(TriangleIteratorIterator a, TriangleIteratorIterator b)
  {
    m_begin = a;
    m_end = b;
    int range = std::distance(m_begin, m_end);

    compute_bbox();

    // sort primitives along longest axis aabb
    sort_primitives();

    if(range > 1)
    {
      // split
      m_left_child  = new AABB_node(m_begin,         m_begin+range/2);
      m_right_child = new AABB_node(m_begin+range/2, m_end);
    }
    else // this is a leaf node, set one primitive
    {
      m_primitive = *m_begin;
    }
  }

  static Triangle_point centroid(TriangleIterator triangle)
  {
    Triangle_point a = (*triangle)[0];
    Triangle_point b = (*triangle)[1];
    Triangle_point c = (*triangle)[2];
    return CGAL::centroid(a,b,c);
  }

  static bool lower_x(const TriangleIterator& f1,
                      const TriangleIterator& f2)
  {
    return centroid(f1).x() < centroid(f2).x();
  }
  static bool lower_y(const TriangleIterator& f1,
                      const TriangleIterator& f2)
  {
    return centroid(f1).y() < centroid(f2).y();
  }
  static bool lower_z(const TriangleIterator& f1,
                      const TriangleIterator& f2)
  {
    return centroid(f1).z() < centroid(f2).z();
  }

  void sort_primitives()
  {
    switch(longest_axis())
    {
      case 0: // sort along x
        std::sort(m_begin,m_end, lower_x);
        break;
      case 1: // sort along y
        std::sort(m_begin,m_end, lower_y);
        break;
      default: // sort along z
        std::sort(m_begin,m_end, lower_z);
    }
  }

  // Return the node's longest axis as 0 (X), 1 (Y) or 2 (Z)
  int longest_axis()
  {
    FT max_size = (std::max)(xsize(),(std::max)(ysize(),zsize()));
    if(max_size == xsize())
      return 0; // axis along x
    if(max_size == ysize())
      return 1; // axis along y
    return 2; // axis along z
  }

  // size of bounding box along each axis
  FT xsize() { return m_bbox.xmax() - m_bbox.xmin(); }
  FT ysize() { return m_bbox.ymax() - m_bbox.ymin(); }
  FT zsize() { return m_bbox.zmax() - m_bbox.zmin(); }

public:

  // --------------------RAY/SEGMENT ORACLES----------------------//

  // Test against primitive if leaf, against box if not
  template<class SegmentType>     ///< tested with Segment_3 and Ray_3
  bool do_intersect(const SegmentType& seg)
  {
    Converter convert;

    if(!is_leaf()) 
    {
      // check against bbox
      return CGAL::do_intersect(convert(seg),m_bbox);
    }
    else 
    {
      // optimization: check first against triangle's bbox
      if (!CGAL::do_intersect(convert(seg),m_bbox))
          return false;

      // check against input primitive 
      return CGAL::do_intersect(Triangle(*m_primitive), seg);
    }
  }

}; // class AABB_node


#endif // AABB_NODE_H_INCLUDED