SPAR/AIModule/BWTA/vendors/CGAL/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h

Go to the documentation of this file.

// Copyright (c) 2006 Tel-Aviv University (Israel).
// 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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h $
// $Id: Arr_polyhedral_sgm.h 50610 2009-07-15 15:21:14Z efif $
// 
// Author(s)     : Efi Fogel         <efif@post.tau.ac.il>

#ifndef CGAL_ARR_POLYHEDRAL_SGM_H
#define CGAL_ARR_POLYHEDRAL_SGM_H

#include <string>
#include <vector>
#include <list>
#include <iostream>

#include <boost/type_traits.hpp>

#include <CGAL/basic.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_traits_with_normals_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/HalfedgeDS_vector.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Aff_transformation_3.h>
#include <CGAL/aff_transformation_tags.h>
#include <CGAL/intersections.h>
#include <CGAL/Polygon_2_algorithms.h>

#include <CGAL/Arr_overlay_2.h>
#include <CGAL/Arr_spherical_gaussian_map_3/Arr_spherical_gaussian_map_3.h>
#include <CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm_traits.h>
#include <CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm_arr_dcel.h>
#include <CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm_overlay.h>
#include <CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm_initializer_visitor.h>

CGAL_BEGIN_NAMESPACE

template <typename PolyhedralSgm,
          typename Polyhedron,
          typename Visitor =
            Arr_polyhedral_sgm_initializer_visitor<PolyhedralSgm, Polyhedron> >
class Arr_polyhedral_sgm_initializer :
  public Arr_sgm_initializer<typename PolyhedralSgm::Base>
{
private:
  // Base type:
  typedef Arr_sgm_initializer<typename PolyhedralSgm::Base> Base;
  typedef typename PolyhedralSgm::Geometry_traits_2         Geometry_traits_2;
  typedef typename Geometry_traits_2::Point_2               Point_2;
  typedef typename Geometry_traits_2::X_monotone_curve_2    X_monotone_curve_2;
  typedef typename Geometry_traits_2::Curve_2               Curve_2;  

  typedef typename Geometry_traits_2::Point_3               Point_3;
  typedef typename Geometry_traits_2::Vector_3              Vector_3;
  
  typedef unsigned int *                                    Coord_index_iter;
  
  // Polyhedron types:
  typedef typename Polyhedron::Vertex_const_handle
    Polyhedron_vertex_const_handle;
  typedef typename Polyhedron::Halfedge_const_handle
    Polyhedron_halfedge_const_handle;

  typedef typename Polyhedron::Vertex_iterator
    Polyhedron_vertex_iterator;
  typedef typename Polyhedron::Halfedge_iterator
    Polyhedron_halfedge_iterator;
  typedef typename Polyhedron::Facet_iterator
    Polyhedron_facet_iterator;

  typedef typename Polyhedron::Halfedge_around_vertex_circulator
    Polyhedron_halfedge_around_vertex_circulator;
  typedef boost::is_same<typename Polyhedron::Plane_3, Vector_3>
    Polyhedron_has_normal;

  struct Normal_equation {
    template <class Facet>
    typename Facet::Plane_3 operator()(Facet & f) {
      typename Facet::Halfedge_handle h = f.halfedge();
      return CGAL::cross_product(h->next()->vertex()->point() -
                                 h->vertex()->point(),
                                 h->next()->next()->vertex()->point() -
                                 h->next()->vertex()->point());
    }
  };

  void compute_planes(Polyhedron & polyhedron, boost::true_type)
  {
    std::transform(polyhedron.facets_begin(), polyhedron.facets_end(),
                   polyhedron.planes_begin(), Normal_equation());
  }
  
  struct Plane_equation {
    template <typename Facet>
    typename Facet::Plane_3 operator()(Facet & f) {
      typename Facet::Halfedge_handle h = f.halfedge();
      return typename Facet::Plane_3(h->vertex()->point(),
                                     h->next()->vertex()->point(),
                                     h->next()->next()->vertex()->point());
    }
  };

  void compute_planes(Polyhedron & polyhedron, boost::false_type)
  {
    std::transform(polyhedron.facets_begin(), polyhedron.facets_end(),
                   polyhedron.planes_begin(), Plane_equation());
  }
  
  template <class HDS, class PointIterator_3>
  class Point_adder {
  private:
    typedef Polyhedron_incremental_builder_3<HDS>               Builder;
    Builder & m_B;

  public:
    typedef typename Polyhedron::Vertex_handle
      Polyhedron_vertex_handle;

    Point_adder(Builder & B) : m_B(B) {}
      
    Polyhedron_vertex_handle operator()(PointIterator_3 pi)
    {
      typedef typename HDS::Vertex      Vertex;
      typedef typename Vertex::Point    Point;
      return m_B.add_vertex(Point((*pi)[0], (*pi)[1], (*pi)[2]));
    }
  };

  template <class HDS> class Point_adder<HDS, Point_3 *> {
  private:
    typedef Polyhedron_incremental_builder_3<HDS>               Builder;
    Builder & m_B;

  public:
    typedef typename Polyhedron::Vertex_handle
      Polyhedron_vertex_handle;

    Point_adder(Builder & B) : m_B(B) {}
      
    Polyhedron_vertex_handle operator()(Point_3 * pi)
    { return m_B.add_vertex(*pi); }
  };
  
  template <class PointIterator_3>
  class Build_surface : public Modifier_base<typename Polyhedron::HalfedgeDS>
  {
  private:
    typedef typename Polyhedron::Vertex_handle
      Polyhedron_vertex_handle;
    typedef typename Polyhedron::Facet_handle
      Polyhedron_facet_handle;
    typedef typename Polyhedron::HalfedgeDS             HDS;
    typedef Polyhedron_incremental_builder_3<HDS>       Builder;
    typedef typename Builder::size_type                 size_type;
    typedef unsigned int *                              Coord_index_iter;

    const PointIterator_3 & m_points_begin;

    const PointIterator_3 & m_points_end;
    
    size_type m_num_points;

    const Coord_index_iter & m_indices_begin;

    const Coord_index_iter & m_indices_end;
    
    size_type m_num_facets;

    unsigned int m_marked_vertex_index;

    unsigned int m_marked_edge_index;

    unsigned int m_marked_facet_index;
    
  public:
    Build_surface(const PointIterator_3 & points_begin,
                  const PointIterator_3 & points_end,
                  unsigned int num_points,
                  const Coord_index_iter & indices_begin,
                  const Coord_index_iter & indices_end,
                  unsigned int num_facets) :
      m_points_begin(points_begin), m_points_end(points_end),
      m_num_points(num_points),
      m_indices_begin(indices_begin), m_indices_end(indices_end),
      m_num_facets(num_facets),
      m_marked_vertex_index(0),
      m_marked_edge_index(0),
      m_marked_facet_index(0)      
    {}

    virtual ~Build_surface() {}

    void set_marked_vertex_index(unsigned int id) {m_marked_vertex_index = id;}

    void set_marked_edge_index(unsigned int id) {m_marked_edge_index = id;}

    void set_marked_facet_index(unsigned int id) {m_marked_facet_index = id;}

    void operator()(HDS & hds)
    {
      // Postcondition: `hds' is a valid polyhedral surface.
      Builder B(hds, true);
      B.begin_surface(m_num_points, m_num_facets);
      // Add the points:
      unsigned int counter = 0;
      Point_adder<HDS, PointIterator_3> add(B);
      for (PointIterator_3 pi = m_points_begin; pi != m_points_end; ++pi) {
        Polyhedron_vertex_handle vh = add(pi);
        if (counter == m_marked_vertex_index) vh->set_marked(true);
        ++counter;
      }
      
      // Add the facets:
      bool facet_ended = true;
      counter = 0;
      for (Coord_index_iter ii = m_indices_begin; ii != m_indices_end; ++ii) {
        int index = *ii;
        if (facet_ended) {
          Polyhedron_facet_handle fh = B.begin_facet();
          if (counter == m_marked_facet_index) fh->set_marked(true);
          B.add_vertex_to_facet(index);
          facet_ended = false;
          continue;
        }
        if (index != -1) {
          B.add_vertex_to_facet(index);
          continue;
        }
        B.end_facet();
        facet_ended = true;
        ++counter;
      }
      B.end_surface();
    }
  };

  Visitor * m_visitor;
  
  unsigned int m_marked_vertex_index;

  unsigned int m_marked_edge_index;

  unsigned int m_marked_facet_index;

  Polyhedron_vertex_const_handle m_src_vertex;

  Polyhedron_vertex_const_handle m_trg_vertex;

  Polyhedron_halfedge_const_handle m_halfedge;

  virtual void handle_new_edge(typename Base::Halfedge_handle edge)
  {
    typedef typename Base::Face_handle          Arr_face_handle;
    typedef typename Base::Vertex_handle        Arr_vertex_handle;
    
    Arr_face_handle src_face = edge->twin()->face();
    Arr_face_handle trg_face = edge->face();
    src_face->set_point(m_src_vertex->point());
    trg_face->set_point(m_trg_vertex->point());

    if (m_visitor) {
      m_visitor->update_dual_vertex(m_src_vertex, src_face);
      m_visitor->update_dual_vertex(m_trg_vertex, trg_face);

      m_visitor->update_dual_halfedge(m_halfedge, edge);
      m_visitor->update_dual_halfedge(m_halfedge, edge->twin());

//       m_visitor->update_dual_face(m_halfedge->opposite()->facet(),
//                                   edge->source());
//       m_visitor->update_dual_face(m_halfedge->facet(), edge->target());
    }
  }
  
  template <class PointIterator_3>  
  void update_polyhedron(Polyhedron & polyhedron,
                         const PointIterator_3 & points_begin,
                         const PointIterator_3 & points_end,
                         unsigned int num_points,
                         const Coord_index_iter indices_begin,
                         Coord_index_iter indices_end,
                         unsigned int num_facets)
  {
    Build_surface<PointIterator_3>
      surface(points_begin, points_end, num_points,
              indices_begin, indices_end, num_facets);
    surface.set_marked_vertex_index(m_marked_vertex_index);
    surface.set_marked_edge_index(m_marked_edge_index);
    surface.set_marked_facet_index(m_marked_facet_index);
    polyhedron.delegate(surface);

    // Mark the marked (half) edges:
    unsigned int counter = 0;
    typedef typename Polyhedron::Edge_iterator Polyhedron_edge_iterator;
    Polyhedron_edge_iterator ei;
    for (ei = polyhedron.edges_begin(); ei != polyhedron.edges_end(); ++ei) {
      if (counter == m_marked_edge_index) {
        // Mark both halfedges:
        ei->set_marked(true);
        ei->opposite()->set_marked(true);
      }
      ++counter;
    }
  }

  template <typename Facet>
  const Vector_3 & get_normal(const Facet & facet, boost::true_type) const
  { return facet->plane(); }
  
  template <typename Facet>
  Vector_3 get_normal(const Facet & facet, boost::false_type) const
  { return facet->plane().orthogonal_vector(); }

  void process_vertex(Polyhedron_vertex_iterator src, bool first_time)
  {
    m_src_vertex = src;

    typedef typename Base::Vertex_handle                Vertex_handle;
    typedef typename Base::Halfedge_handle              Halfedge_handle;

    Vertex_handle invalid_vertex;
    
    // For each vertex, traverse incident faces:
    Polyhedron_halfedge_around_vertex_circulator hec = src->vertex_begin();

    // If the vertex is not a real vertex of the polyhedron, advance to the
    // next one:
    if (circulator_size(hec) == 0) {
      process_vertex(++src, first_time);
      return;
    }

    CGAL_assertion(circulator_size(hec) >= 3);
    Polyhedron_halfedge_around_vertex_circulator begin_hec = hec;
    Polyhedron_halfedge_around_vertex_circulator next_hec = hec;
    ++next_hec;

    /* If this is not the first invocation, advance the halfedge iterator
     * until its source vertex is processed. It is guaranteed to reach such
     * a halfedge on consecutive invocations.
     */
    if (!first_time) {
      while (!hec->opposite()->vertex()->processed()) {
        hec = next_hec;
        begin_hec = hec;
        ++next_hec;
      }
    }

    // Traverse the incident halfedges:
    do {
      if (!next_hec->processed()) {

        Vector_3 normal1 =
          get_normal(hec->facet(), Polyhedron_has_normal());
        Vector_3 normal2 =
          get_normal(next_hec->facet(), Polyhedron_has_normal());

        m_trg_vertex = next_hec->opposite()->vertex();

#if 0
        std::cout << "process_vertex trg: "
                  << static_cast<float>(todouble(m_trg_vertex->point().x()))
                  << ","
                  << static_cast<float>(todouble(m_trg_vertex->point().y()))
                  << ","
                  << static_cast<float>(todouble(m_trg_vertex->point().z()))
                  << std::endl;
#endif
        
        m_halfedge = next_hec;
#if 0
        Halfedge_handle he = this->insert(normal1, normal2);
#else
        Vertex_handle v1 = hec->facet()->vertex();
        Vertex_handle v2 = next_hec->facet()->vertex();
        /* The arc might be non-x-monotone. In this case, it is broken into 2
         * x-monotone curves. The halfedges of both are obtained.
         */
        typedef typename std::list<Halfedge_handle>     Halfedge_list;
        typedef typename Halfedge_list::iterator        Halfedge_list_iter;
        Halfedge_list hes;
        if (first_time) {
          this->insert(normal1, normal2, std::back_inserter(hes));
          Halfedge_list_iter first = hes.begin();
          Halfedge_list_iter last = hes.end();
          --last;
          hec->facet()->set_vertex((*first)->source());
          next_hec->facet()->set_vertex((*last)->target());
          first_time = false;
        } else {
          if (v1 != invalid_vertex && v2 != invalid_vertex) {
            this->insert(normal1, v1, normal2, v2, std::back_inserter(hes));
          } else if (v1 != invalid_vertex) {
            this->insert(normal1, v1, normal2, std::back_inserter(hes));
            Halfedge_list_iter last = hes.end();
            --last;
            next_hec->facet()->set_vertex((*last)->target());
          } else if (v2 != invalid_vertex) {
            this->insert(normal1, normal2, v2, std::back_inserter(hes));
            Halfedge_list_iter first = hes.begin();
            hec->facet()->set_vertex((*first)->source());
          } else CGAL_error();
        }
#endif
        next_hec->set_processed(true);
        next_hec->opposite()->set_processed(true);
        Halfedge_list_iter first = hes.begin();
        if (v1 != invalid_vertex && v2 != invalid_vertex)
          handle_new_edge(*first);

        if (m_visitor) {
          m_visitor->update_dual_face(m_halfedge->opposite()->facet(),
                                      (*first)->source());
          // m_visitor->update_dual_face(m_halfedge->facet(), (*first)->target());
        }
      }
      hec = next_hec;
      ++next_hec;
    } while (hec != begin_hec);
    src->set_processed(true);

    // Traverse recursively:
    hec = src->vertex_begin();
    begin_hec = hec;
    do {
      if (!(hec->opposite()->vertex()->processed()))
        process_vertex(hec->opposite()->vertex(), false);
      ++hec;
    } while (hec != begin_hec);
  }
  
  void compute_sgm(Polyhedron & polyhedron)
  {
    typedef typename Base::Vertex_handle                Vertex_handle;

    // Clear the polyhedron:
    Polyhedron_facet_iterator fi;
    for (fi = polyhedron.facets_begin(); fi != polyhedron.facets_end(); ++fi)
      fi->set_vertex(Vertex_handle());

    Polyhedron_halfedge_iterator hei;
    for (hei = polyhedron.halfedges_begin(); hei != polyhedron.halfedges_end();
         ++hei)
      hei->set_processed(false);

    Polyhedron_vertex_iterator vi;
    for (vi = polyhedron.vertices_begin(); vi != polyhedron.vertices_end();
         ++vi)
      vi->set_processed(false);

    // Traverse all verticess recursively:
    process_vertex(polyhedron.vertices_begin(), true);
  }
  
public:
  Arr_polyhedral_sgm_initializer(PolyhedralSgm & sgm) :
    Base(sgm),
    m_visitor(NULL),
    m_marked_vertex_index(0),
    m_marked_edge_index(0),
    m_marked_facet_index(0)
  {}
    
  virtual ~Arr_polyhedral_sgm_initializer() {}

  void operator()(Polyhedron & polyhedron, Visitor * visitor = NULL)
  {
#if 0
    std::copy(polyhedron.points_begin(), polyhedron.points_end(),
              std::ostream_iterator<Point_3>(std::cout, "\n"));
#endif

    m_visitor = visitor;

#if 0
    if (!polyhedron.normalized_border_is_valid())
      polyhedron.normalize_border();
#else
    polyhedron.normalize_border();
#endif
#if 1
    compute_planes(polyhedron, Polyhedron_has_normal());
#endif

    compute_sgm(polyhedron);
  }
  
  template <class PointIterator_3>
  void operator()(const PointIterator_3 & points_begin,
                  const PointIterator_3 & points_end,
                  unsigned int num_points,
                  const Coord_index_iter indices_begin,
                  Coord_index_iter indices_end,
                  unsigned int num_facets,
                  Visitor * visitor = NULL)
  {
    m_visitor = visitor;
 
    Polyhedron polyhedron;
    update_polyhedron(polyhedron, points_begin, points_end, num_points,
                      indices_begin, indices_end, num_facets);

#if 0
    std::copy(polyhedron.points_begin(), polyhedron.points_end(),
              std::ostream_iterator<Point_3>(std::cout, "\n"));
#endif

#if 0
    if (!polyhedron.normalized_border_is_valid())
      polyhedron.normalize_border();
#else
    polyhedron.normalize_border();
#endif
#if 1
    compute_planes(polyhedron, Polyhedron_has_normal());
#endif

    compute_sgm(polyhedron);
    polyhedron.clear();
  }

  void set_marked_vertex_index(unsigned int id) {m_marked_vertex_index = id;}

  void set_marked_edge_index(unsigned int id) {m_marked_edge_index = id;}

  void set_marked_facet_index(unsigned int id) {m_marked_facet_index = id;}
};

template <class Geometry_traits_2,
#ifndef CGAL_CFG_NO_TMPL_IN_TMPL_PARAM
          template <class T>
#endif
          class T_Dcel = Arr_polyhedral_sgm_arr_dcel>
class Arr_polyhedral_sgm :
  public Arr_spherical_gaussian_map_3<Geometry_traits_2, T_Dcel>
{
private:
  typedef Arr_polyhedral_sgm<Geometry_traits_2, T_Dcel>     Self;
  
public:
  typedef typename Geometry_traits_2::Point_3               Point_3;
  typedef typename Geometry_traits_2::Vector_3              Vector_3;

#ifndef CGAL_CFG_NO_TMPL_IN_TMPL_PARAM
  typedef T_Dcel<Geometry_traits_2>                         Dcel;
#else
  typedef typename T_Dcel::template Dcel<Geometry_traits_2> Dcel;
#endif
  
  // For some reason MSVC barfs on the friend statement below. Therefore,
  // we declare the Base to be public to overcome the problem.
  typedef Arr_spherical_gaussian_map_3<Geometry_traits_2, T_Dcel>   Base;

  typedef Arr_polyhedral_sgm_overlay<Self>
    Arr_polyhedral_sgm_overlay;

#if 0

  template <class Polyhedron, class Visitor>
  friend class Arr_polyhedral_sgm_initializer<Self, Polyhedron, Visitor>;
#endif
  
private:
  Point_3 m_center;
  
  bool m_dirty_center;
  
  void calculate_center()
  {
    // Count them:
    typename Base::Face_handle fi;
    for (fi = this->faces_begin(); fi != this->faces_end(); fi++) {
      const Point_3 & p = fi->point();
      Vector_3 v = p - CGAL::ORIGIN;
      m_center = m_center + v;
    }

    typename Base::Size num = this->number_of_faces();
    m_center =
      Point_3(m_center.x() / num, m_center.y() / num, m_center.z() / num);

    m_dirty_center = false;
  }
  
public:
  Arr_polyhedral_sgm() : m_dirty_center(true) {}
  
  Arr_polyhedral_sgm(const Self & sgm)
  {
    assign(sgm);
  }

  void assign(const Self & sgm)
  {
    // Call the assign of the base class.
    Base::assign(sgm);

    typename Dcel::Face_iterator fit;
    typename Base::Face_const_iterator fit1 = sgm.faces_begin();
          
    // Set the points of the faces.
    for (fit = this->_dcel().faces_begin(); fit != this->_dcel().faces_end();
         ++fit)
    {
      fit->set_point (fit1->point());
      ++fit1;
    }
                  
    typename Dcel::Edge_iterator eit;
    typename Base::Edge_const_iterator eit1 = sgm.edges_begin();

    // Set the arr_mask of the edges
    for (eit = this->_dcel().edges_begin(); eit != this->_dcel().edges_end();
         ++eit)
    {
      eit->set_arr(eit1->arr_mask());
      eit->opposite()->set_arr(eit1->arr_mask());
      ++eit1;
    }
  }
  
  virtual ~Arr_polyhedral_sgm() { clear(); }

  void clear()
  {
    m_dirty_center = true;
    Base::clear();
  }
  
 // /*! Compute the minkowski sum of a range of objects of type
 //  * Arr_polyhedral_sgm
 //  */
 // template <class SgmIterator>  
 // void minkowski_sum(SgmIterator begin, SgmIterator end)
 // {
        //typename SgmIterator::value_type * sgm1 = *begin++;
 //   typename SgmIterator::value_type * sgm2 = *begin;
 //   minkowski_sum(sgm1, sgm2);
 // }

 // /*! Compute the minkowski sum of a range of objects of type
 //  * Arr_polyhedral_sgm
 //  */
 // template <typename SgmIterator, typename OverlayTraits>  
 // void minkowski_sum(SgmIterator begin, SgmIterator end,
 //                    OverlayTraits & overlay_traits)
 // {
 //   typename SgmIterator::value_type * sgm1 = *begin++;
 //   typename SgmIterator::value_type * sgm2 = *begin;
 //   minkowski_sum(sgm1, sgm2, overlay_traits);
 // }
 
  template <class Arr_polyhedral_sgm>
  void minkowski_sum(const Arr_polyhedral_sgm & sgm1,
                     const Arr_polyhedral_sgm & sgm2)
  {
    // Compute the overlays:
    Arr_polyhedral_sgm_overlay sgm_overlay;
    CGAL::overlay(sgm1, sgm2, *this, sgm_overlay);
    // print_stat();
  }

  template <class Arr_polyhedral_sgm, typename OverlayTraits>
  void minkowski_sum(const Arr_polyhedral_sgm & sgm1,
                     const Arr_polyhedral_sgm & sgm2,
                     OverlayTraits & overlay_traits)
  { CGAL::overlay(sgm1, sgm2, *this, overlay_traits); }
  
  unsigned int number_of_vertices() const
  { return (static_cast<const Base*>(this))->number_of_faces(); }
  
  unsigned int number_of_edges() const
  {
    unsigned int size = 0;
    typename Base::Vertex_const_iterator vit;
    for (vit = this->vertices_begin(); vit != this->vertices_end(); ++vit)
      if (vit->degree() == 2) size++;
    return (static_cast<const Base*>(this))->number_of_edges() - size;
  }

  unsigned int number_of_facets() const
  {
    unsigned int size = 0;
    typename Base::Vertex_const_iterator vit;
    for (vit = this->vertices_begin(); vit != this->vertices_end(); ++vit)
      if (vit->degree() > 2) size++;
    return size;
  }
  
  void print_stat()
  {
    Base::print_stat();
    
    std::cout << "Polyhedron"
              << ", no. facets: " << number_of_facets()
              << ", no. edges: " << number_of_edges()
              << ", no. vertices: " << number_of_vertices()
              << std::endl;
  }
};

CGAL_END_NAMESPACE

#endif