CGAL users discussion list

["Laurent Rineau (GeometryFactory)" <>]

Le mercredi 08 février 2012 09:05:40 Arnaud A. a écrit :
> Please find attached four files: two polyhedral files (conteneur.off,
> copeau.off) and two 3D mesh files (conteneur.mesh,copeau.mesh). I would like
> to make a cut (boolean operation) between conteneur.mesh and copeau.mesh,
> and remesh these two domains in order to have a compatible mesh at the
> interface.
> 
> Is it possible ?

You said the boundaries of your two polyhedra intersect, but they do not. In 
that case that is more simple.


I attach a file "Multi_domain_*.h", that defines a new class, a variation of 
CGAL::Polyhedral_mesh_domain. It is a prototype class, that is not yet in 
CGAL, and will not be either in the next release of CGAL.

Use it that way:

// The typedefs
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::Type Polyhedron;
typedef CGAL::Multi_domain_polyhedral_mesh_domain_3<Polyhedron, K> 
Mesh_domain;

[...]

// Fill the domain

Mesh_domain mesh_domain; // empty domain

Polyhedron conteneur; // empty polyhedron
std::ifstream in("conteneur.off");
in >> conteneur; // read the first polyhedron
if(!in) return EXIT_FAILURE;
in.close();

Polyhedron copeau; // second polyhedron, empty at first
in.open("copeau.off");
in >> copeau; // read the second polyhedron
if(!in) return EXIT_FAILURE;
in.close();

mesh_domain.add_subdomain(conteneur, 1, 60);
// index == 1 for that domain, and detect sharp edges for angle>60°

mesh_domain.add_subdomain(copeau, 1);
// index == 2, not sharp edges

mesh_domain.build(); // finish the creation of the domain object



Then use 'mesh_domain' as the usual mesh domain variable. Note that it may 
require CGAL-3.9.

-- 
Laurent Rineau, PhD
R&D Engineer at GeometryFactory           http://www.geometryfactory.com/
Release Manager of the CGAL Project       http://www.cgal.org/
// Copyright (c) 2011 GeometryFactory Sarl (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: https://geometryfactory.svn.codebasehq.com/hutchinson/code/mesh_generation/Multi_domain_polyhedral_mesh_domain_3.h $
// $Id: Multi_domain_polyhedral_mesh_domain_3.h 24 2011-12-16 15:42:21Z lrineau $
//
// Author(s)     : Laurent Rineau

#ifndef CGAL_MULTI_DOMAIN_POLYHEDRAL_MESH_DOMAIN_3_H
#define CGAL_MULTI_DOMAIN_POLYHEDRAL_MESH_DOMAIN_3_H

#include <CGAL/Polyhedral_mesh_domain_3.h>

#include <CGAL/Mesh_domain_with_polyline_features_3.h>

#include <CGAL/Mesh_3/Detect_polylines_in_polyhedra.h>
#include <CGAL/Mesh_3/Polyline_with_context.h>
#include <CGAL/Mesh_3/Detect_features_in_polyhedra.h>

#include "angle_3_extension.h"

#include <boost/foreach.hpp>

namespace CGAL {

namespace Mesh_3 {

template<typename TriangleAccessor, 
         typename Gt, 
         typename Subdomain_index = int>
class Multi_domain_polyhedral_primitive
{
public:
  typedef Multi_domain_polyhedral_primitive<
    TriangleAccessor,
    Gt,
    Subdomain_index>                                  Self;
  typedef const Self*                                 Id;
  typedef typename TriangleAccessor::Triangle_handle  Handle;
  typedef typename Gt::Triangle_3                     Datum;
  typedef typename Gt::Point_3                        Point;

  Multi_domain_polyhedral_primitive(const Handle& h, 
                                    const Subdomain_index subdomain_index = 1)
    : handle_(h)
    , subdomain_index_(subdomain_index)
  {}

  Datum datum() const
  {
    return TriangleAccessor().triangle(handle_);
  }

  Id id() const
  {
    return this;
  }

  Handle handle() const {
    return handle_;
  }

  Subdomain_index subdomain_index() const {
    return subdomain_index_;
  }

  Point reference_point() const
  {
    return TriangleAccessor().triangle(handle_).vertex(0);
  }

private:
  Handle handle_;
  Subdomain_index subdomain_index_;
}; // end class template Multi_domain_polyhedral_primitive

template <typename Gt, typename Angle_3>
class First_is_closest_along_ray
{
  typedef typename Gt::Ray_3     Ray;
  typedef typename Gt::Vector_3  Vector;
  typedef typename Gt::Point_3   Point;

  const Vector v;

public:
  First_is_closest_along_ray(const Vector& v)
    : v(v)
  {}

  template <typename T>
  bool operator()(const T& a, const T& b) const 
  {
    return Angle_3()(a.first, b.first, v) == ACUTE;
  }
}; // end First_is_closest_along_ray

} // end namespace Mesh_3

template<typename Polyhedron_,
         typename IGT_,
         typename Subdomain_index_ = int,
         typename TriangleAccessor=Triangle_accessor_3<Polyhedron_,IGT_>
         >
class Multi_domain_polyhedral_mesh_domain_base_3 
{
public:
  typedef CGAL::Tag_true Use_exact_intersection_construction_tag;
  typedef typename Mesh_3::details::IGT_generator<
    IGT_,Use_exact_intersection_construction_tag>::type IGT;

  typedef IGT R;

  /// Geometric object types
  typedef typename IGT::Point_3    Point_3;
  typedef typename IGT::Segment_3  Segment_3;
  typedef typename IGT::Ray_3      Ray_3;
  typedef typename IGT::Line_3     Line_3;

  typedef typename IGT::FT         FT;

  typedef Polyhedron_ Polyhedron;

  //-------------------------------------------------------
  // Index Types
  //-------------------------------------------------------
  /// Type of indexes for cells of the input complex
  typedef Subdomain_index_                                Subdomain_index;
  typedef boost::optional<Subdomain_index>                Subdomain;
  /// Type of indexes for surface patch of the input complex
  typedef Subdomain_index                                 Surface_patch_index;
  typedef boost::optional<Surface_patch_index>            Surface_patch;

  typedef int                                             Curve_segment_index;
  typedef int                                             Corner_index;
  typedef CGAL::Tag_true                                  Has_features;

  /// Type of indexes to characterize the lowest dimensional face of the input
  /// complex on which a vertex lie
  typedef Subdomain_index                                 Index;

  typedef CGAL::cpp0x::tuple<Point_3,Index,int> Intersection;


private:
  typedef Mesh_3::Multi_domain_polyhedral_primitive<
  TriangleAccessor, IGT, Subdomain_index>                 AABB_primitive;
  typedef class AABB_traits<IGT,AABB_primitive>           AABB_traits;
  typedef class CGAL::AABB_tree<AABB_traits>              AABB_tree;
private:
  typedef typename AABB_tree::Primitive_id                AABB_primitive_id;
  typedef typename AABB_tree::Primitive                   Primitive;
  typedef typename AABB_traits::Bounding_box              Bounding_box;
  
private:
  /// The AABB tree: intersection detection and more
  AABB_tree tree_;
  std::vector<std::pair<Point_3, int> > initial_points;

  /// Allowed query types
  typedef boost::mpl::vector<Segment_3, Ray_3, Line_3> Allowed_query_types;

  // cache queries and intersected primitive
  typedef typename boost::make_variant_over<Allowed_query_types>::type Cached_query;
  mutable bool has_cache;
  mutable Cached_query cached_query;
  mutable AABB_primitive_id cached_primitive_id;

public:
  /// Default constructor
  Multi_domain_polyhedral_mesh_domain_base_3()
    : tree_()
    , has_cache(false) {}

  void add_subdomain(Polyhedron& p, int index// , 
                     // double detect_sharp_features_angle = 0.
                     ) 
  {
    // if(detect_sharp_features_angle != 0.) {
    //   detect_features(p, detect_sharp_features_angle);
    // }
    std::set<Point_3> all_points_set;
    for(typename TriangleAccessor::Triangle_iterator 
          it = TriangleAccessor().triangles_begin(p),
          end = TriangleAccessor().triangles_end(p);
        it != end; ++it)
    {
      all_points_set.insert(TriangleAccessor().triangle(it).vertex(0));
      all_points_set.insert(TriangleAccessor().triangle(it).vertex(1));
      all_points_set.insert(TriangleAccessor().triangle(it).vertex(2));
      tree_.insert(Primitive(it, index));
    }
    std::vector<Point_3> all_points(all_points_set.begin(),
                                    all_points_set.end());

    const typename Polyhedron::size_type initial_points_nb = 
      (std::min)(all_points.size(), 8LU);

    std::random_shuffle(all_points.begin(), all_points.end());
    for(typename Polyhedron::size_type i = 0;
        i <= initial_points_nb; ++i)
    {
      initial_points.push_back(std::make_pair(all_points[i], index));
    }
  }

  void build() {
    tree_.build();
  }
  
  /// Destructor
  ~Multi_domain_polyhedral_mesh_domain_base_3() { 
  }

  /**
   * Constructs  a set of \ccc{n} points on the surface, and output them to
   *  the output iterator \ccc{pts} whose value type is required to be
   *  \ccc{std::pair<Points_3, Index>}.
   */
  struct Construct_initial_points
  {
    Construct_initial_points(const Multi_domain_polyhedral_mesh_domain_base_3& domain)
      : r_domain_(domain) {}

    template<class OutputIterator>
    OutputIterator operator()(OutputIterator pts) const {
      return std::copy(r_domain_.initial_points.begin(), 
                       r_domain_.initial_points.end(), 
                       pts);
    }

  private:
    const Multi_domain_polyhedral_mesh_domain_base_3& r_domain_;
  };

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


  /**
   * Returns true if point~\ccc{p} is in the domain. If \ccc{p} is in the
   *  domain, the parameter index is set to the index of the subdomain
   *  including $p$. It is set to the default value otherwise.
   */
  struct Is_in_domain
  {
    Is_in_domain(const Multi_domain_polyhedral_mesh_domain_base_3& domain)
      : r_domain_(domain) {}

    Subdomain operator()(const Point_3& p) const;
  private:
    const Multi_domain_polyhedral_mesh_domain_base_3& r_domain_;
  };

  Is_in_domain is_in_domain_object() const { return Is_in_domain(*this); }

  Point_3 project_on_surface(const Point_3& p) const
  {
    return tree_.closest_point(p);
  }
  
  /**
   * Returns true is the element \ccc{type} intersect properly any of the
   * surface patches describing the either the domain boundary or some
   * subdomain boundary.
   * \ccc{Type} is either \ccc{Segment_3}, \ccc{Ray_3} or \ccc{Line_3}.
   * Parameter index is set to the index of the intersected surface patch
   * if \ccc{true} is returned and to the default \ccc{Surface_patch_index}
   * value otherwise.
   */
  struct Do_intersect_surface
  {
    Do_intersect_surface(const Multi_domain_polyhedral_mesh_domain_base_3& domain)
      : r_domain_(domain) {}

    template <typename Query>
    typename boost::enable_if<typename boost::mpl::contains<Allowed_query_types,
                                                            Query>::type,
                              Surface_patch>::type
    operator()(const Query& q) const
    {
      boost::optional<AABB_primitive_id> primitive_id = r_domain_.tree_.any_intersected_primitive(q);
      if ( primitive_id )
      { 
        r_domain_.cache_primitive(q, *primitive_id);
        return Surface_patch(r_domain_.make_surface_index(*primitive_id));
      } else {      
        return Surface_patch();
      }
    }

  private:
    const Multi_domain_polyhedral_mesh_domain_base_3& r_domain_;
  };

  Do_intersect_surface do_intersect_surface_object() const
  {
    return Do_intersect_surface(*this);
  }

  /**
   * Returns a point in the intersection of the primitive \ccc{type}
   * with some boundary surface.
   * \ccc{Type1} is either \ccc{Segment_3}, \ccc{Ray_3} or \ccc{Line_3}.
   * The integer \ccc{dimension} is set to the dimension of the lowest
   * dimensional face in the input complex containing the returned point, and
   * \ccc{index} is set to the index to be stored at a mesh vertex lying
   * on this face.
   */
  struct Construct_intersection
  {
    Construct_intersection(const Multi_domain_polyhedral_mesh_domain_base_3& domain)
      : r_domain_(domain) {}

    template <typename Query>
    typename boost::enable_if<typename boost::mpl::contains<Allowed_query_types,
                                                            Query>::type,
                              Intersection>::type
    operator()(const Query& q) const
    {
      typedef boost::optional<
        typename AABB_tree::Object_and_primitive_id> AABB_intersection;
      typedef Point_3 Bare_point;

      AABB_intersection intersection;
      if(r_domain_.query_is_cached(q))
      {
        const AABB_primitive_id primitive_id = r_domain_.cached_primitive_id;
        Object o = IGT().intersect_3_object()(primitive_id->datum(),
                                              q);
        intersection = AABB_intersection(std::make_pair(o, primitive_id));
      } else {

      CGAL_precondition(r_domain_.do_intersect_surface_object()(q));

      intersection = r_domain_.tree_.any_intersection(q);
      }
      if ( intersection )
      {
        // Get primitive
        AABB_primitive_id primitive_id = (*intersection).second;
        
        // intersection may be either a point or a segment
        if ( const Bare_point* p_intersect_pt =
                              object_cast<Bare_point>(&(*intersection).first) )
        {
          return Intersection(*p_intersect_pt,
                              r_domain_.index_from_surface_patch_index(
                                r_domain_.make_surface_index(primitive_id)),
                              2);
        }
        else if ( const Segment_3* p_intersect_seg =
                              object_cast<Segment_3>(&(*intersection).first) )
        {
          return Intersection(p_intersect_seg->source(),
                              r_domain_.index_from_surface_patch_index(
                                r_domain_.make_surface_index(primitive_id)),
                              2);
        }
        else
          CGAL_error_msg("Mesh_3 error : AABB_tree any_intersection result is "
                         "not a point nor a segment");
      }

      // Should not happen
      return Intersection();
    }

  private:
    const Multi_domain_polyhedral_mesh_domain_base_3& r_domain_;
  };

  Construct_intersection construct_intersection_object() const
  {
    return Construct_intersection(*this);
  }
  
  
  /**
   * Returns the index to be stored in a vertex lying on the surface identified
   * by \c index.
   */
  Index index_from_surface_patch_index(const Surface_patch_index& index) const
  { return Index(index); }

  /**
   * Returns the index to be stored in a vertex lying in the subdomain
   * identified by \c index.
   */
  Index index_from_subdomain_index(const Subdomain_index& index) const
  { return Index(index); }

  /**
   * Returns the \c Surface_patch_index of the surface patch
   * where lies a vertex with dimension 2 and index \c index.
   */
  Surface_patch_index surface_patch_index(const Index& index) const
  { return index; }

  /**
   * Returns the index of the subdomain containing a vertex
   *  with dimension 3 and index \c index.
   */
  Subdomain_index subdomain_index(const Index& index) const
  { return index; }
  
public:
  Surface_patch_index make_surface_index(
    const AABB_primitive_id& primitive_id = AABB_primitive_id() ) const
  {
    return primitive_id->subdomain_index();
  }

  // Undocumented function, used to implement a sizing field that
  // computes lfs using this AABB tree. That avoids to rebuild the same
  // tree.
  const AABB_tree& aabb_tree() const {
    return tree_;
  }

protected:
  void add_primitives(const Polyhedron& p)
  {
    tree_.insert(TriangleAccessor().triangles_begin(p),
                 TriangleAccessor().triangles_end(p));
    
    tree_.build();
  }

private:
  template <typename Query>
  void cache_primitive(const Query& q, 
                       const AABB_primitive_id id) const
  {
    cached_query = Cached_query(q);
    has_cache = true;
    cached_primitive_id = id;
  }

  template <typename Query>
  bool query_is_cached(const Query& q) const {
    return has_cache && (cached_query == Cached_query(q));
  }

private:
  // Disabled copy constructor & assignment operator
  typedef Multi_domain_polyhedral_mesh_domain_base_3 Self;
  Multi_domain_polyhedral_mesh_domain_base_3(const Self& src);
  Self& operator=(const Self& src);

};  // end class Multi_domain_polyhedral_mesh_domain_base_3

template<typename P_, typename IGT_, typename TA, typename Tag>
typename Multi_domain_polyhedral_mesh_domain_base_3<P_,IGT_,TA,Tag>::Subdomain
Multi_domain_polyhedral_mesh_domain_base_3<P_,IGT_,TA,Tag>::
Is_in_domain::operator()(const Point_3& p) const
{
  // Shoot ray
  typename IGT::Construct_ray_3 ray = IGT().construct_ray_3_object();
  typename IGT::Construct_vector_3 vector = IGT().construct_vector_3_object();
  
  Random_points_on_sphere_3<Point_3> random_point(1.);

  const typename IGT::Vector_3 v = vector(CGAL::ORIGIN,*random_point);

  const Ray_3 ray_shot = ray(p, v);

  typedef typename AABB_tree::Object_and_primitive_id Object_and_primitive_id;

  std::vector<Object_and_primitive_id> all_intersections;
  r_domain_.tree_.all_intersections(ray_shot, 
                                    std::back_inserter(all_intersections));
  if(all_intersections.empty()) return Subdomain();

  typedef Mesh_3::First_is_closest_along_ray<IGT, 
                Filtered_angle_3_ext<IGT> > First_along_v;
  First_along_v first_along_v(v);

  typedef std::pair<Point_3, int> Pair_point_and_index;
  std::set< Pair_point_and_index,
            First_along_v > intersections(first_along_v);

  // std::cerr << "\nintersections: ";
  BOOST_FOREACH(const Object_and_primitive_id& obj_pri_id, 
                all_intersections)
  {
    if(const Point_3* p = object_cast<Point_3>(&obj_pri_id.first))
    {
      // std::cerr << obj_pri_id.second->subdomain_index() << " ";
      intersections.insert(std::make_pair(*p, 
                                          obj_pri_id.second->subdomain_index()));
    }
  }
  // std::cerr << "\n";

  std::stack<int> indices_stack;

  BOOST_REVERSE_FOREACH( const Pair_point_and_index& p_and_index,
                         intersections )
  {
    const int& index = p_and_index.second;
    if( !indices_stack.empty() && indices_stack.top() == index )
      indices_stack.pop();
    else 
      indices_stack.push(index);
  }

  if(indices_stack.empty()) {
    // std::cerr << "result empty\n";
    return Subdomain();
  } else {
    // std::cerr << "result: " << indices_stack.top() << std::endl;
    return Subdomain(Subdomain_index(indices_stack.top()));
  }
}

template<typename Polyhedron_,
         typename IGT_,
         typename Subdomain_index_ = int,
         typename TriangleAccessor=Triangle_accessor_3<Polyhedron_,IGT_>
         >
class Multi_domain_polyhedral_mesh_domain_3 :
  public Mesh_domain_with_polyline_features_3<
    Multi_domain_polyhedral_mesh_domain_base_3<Polyhedron_,
                                               IGT_,
                                               Subdomain_index_,
                                               TriangleAccessor> 
  >
{
  typedef Multi_domain_polyhedral_mesh_domain_base_3<Polyhedron_,
                                                     IGT_,
                                                     Subdomain_index_,
                                                     TriangleAccessor> BaseBase;
  typedef Mesh_domain_with_polyline_features_3<BaseBase> Base;

public:
  typedef Polyhedron_ Polyhedron;
  typedef typename IGT_::FT FT;
  typedef typename IGT_::Point_3 Point_3;
  typedef typename Base::Surface_patch_index Surface_patch_index;
  typedef typename Base::Curve_segment_index Curve_segment_index;

  /// Default constructor
  Multi_domain_polyhedral_mesh_domain_3()
    : Base(), maximal_index(0) {}

  /// Detect features
  void detect_features(Polyhedron& polyhedron_, FT angle_in_degree)
  {
    // Get sharp features
    Mesh_3::detect_features(polyhedron_,angle_in_degree);
  
    // Get polylines
    typedef std::vector<Point_3> Bare_polyline;
    typedef Mesh_3::Polyline_with_context<Surface_patch_index, Curve_segment_index,
                                          Bare_polyline > Polyline;
  
    std::vector<Polyline> polylines;
    typedef std::back_insert_iterator<std::vector<Polyline> > Output_iterator;

    Mesh_3::detect_polylines<Polyhedron,Polyline,Output_iterator>(
                                                                  &polyhedron_, std::back_inserter(polylines));
    
    // Insert polylines in domain
    Mesh_3::Extract_bare_polyline<Polyline> extractor;
  
    this->add_features(
      boost::make_transform_iterator(polylines.begin(),extractor),
      boost::make_transform_iterator(polylines.end(),extractor));
  }
  
  void add_subdomain(Polyhedron& p, int index, 
                     double detect_sharp_features_angle = 0.) 
  {
    if(detect_sharp_features_angle != 0.) {
      detect_features(p, detect_sharp_features_angle);
    }
    Base::add_subdomain(p, index);
    if(maximal_index < index) maximal_index = index;
  }

  struct Construct_initial_points
  {
    Construct_initial_points(const Multi_domain_polyhedral_mesh_domain_3& domain)
      : r_domain_(domain) {}

    template<class OutputIterator>
    OutputIterator operator()(OutputIterator pts) const {
      return
        this->r_domain_.Base::construct_initial_points_object()(pts);
    }

  private:
    const Multi_domain_polyhedral_mesh_domain_3& r_domain_;
  };

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

  int maximal_index;

}; // end class Multi_domain_polyhedral_mesh_domain_with_features_3

} // end namespace CGAL

#endif // CGAL_MULTI_DOMAIN_POLYHEDRAL_MESH_DOMAIN_3_H