surface_mesher.h

很多二维 三维几何计算算法 C++ 类库

// Copyright (c) 2003-2007  INRIA Sophia-Antipolis (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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Surface_mesher/include/CGAL/Surface_mesher/Surface_mesher.h $
// $Id: Surface_mesher.h 38107 2007-04-12 15:50:19Z lrineau $
//
//
// Author(s) : Steve Oudot,
//             David Rey,
//             Mariette Yvinec,
//             Laurent Rineau,
//             Andreas Fabri



#ifndef CGAL_SURFACE_MESHER_SURFACE_MESHER_H
#define CGAL_SURFACE_MESHER_SURFACE_MESHER_H

#include <CGAL/Mesher_level.h>
#include <CGAL/Mesh_2/Triangulation_mesher_level_traits_3.h>
#include <CGAL/Double_map.h>
#include <list>
#include <string>
#include <sstream>

namespace CGAL {

  namespace Surface_mesher {

  // NB: by convention, the priority queue is sorted with respect to the
  // first element of the list of criteria

  template <
    class C2T3,
    class Surface_,
    class SurfaceMeshTraits,
    class Criteria_
    >
  class Surface_mesher_base
    : public Triangulation_mesher_level_traits_3<typename C2T3::Triangulation>
  {
  public:
    typedef C2T3 Complex_2_in_triangulation_3;
    typedef Surface_ Surface;
    typedef SurfaceMeshTraits Surface_mesh_traits;
    typedef Criteria_ Criteria;

    typedef typename C2T3::Triangulation Tr;
    typedef typename Tr::Point Point;
    typedef typename Tr::Edge Edge;
    typedef typename Tr::Vertex_handle Vertex_handle;
    typedef typename Tr::Cell_handle Cell_handle;

    typedef typename Tr::Geom_traits GT;

    typedef Triangulation_mesher_level_traits_3<Tr> Triangulation_traits;
    typedef typename Triangulation_traits::Zone Zone;

    typedef typename Tr::Finite_facets_iterator Finite_facets_iterator;
    typedef typename Tr::Facet Facet;

    typedef typename Criteria::Quality Quality;

    typedef Double_map<Facet, Quality> Bad_facets;

    // Constructor
    Surface_mesher_base (C2T3& co, 
                         const Surface& s, 
                         const Surface_mesh_traits& mesh_traits,
                         const Criteria& c) :
      Triangulation_mesher_level_traits_3<Tr>(co.triangulation()),
      c2t3(co),
      tr(co.triangulation()),
      surf(s),
      meshtraits(mesh_traits),
      criteria(c)
    {
#ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
      std::cerr << "CONS: Surface_mesher_base\n";
#endif
    }

  protected:
    C2T3& c2t3;
    Tr& tr;     // Associated triangulation reference
    const Surface& surf;  // Surface
    const Surface_mesh_traits& meshtraits; // Surface mesh traits
    const Criteria& criteria;  // Meshing criteria
    Bad_facets facets_to_refine;  // Set of facets to refine

  public:

    // Helper functions
    Facet mirror_facet(const Facet& f) const
    {
      return tr.mirror_facet(f);
    }

    static void set_facet_visited(Facet f)
    {
      f.first->set_facet_visited(f.second);
    }

    static void set_facet_surface_center(Facet f, const Point& p)
    {
      f.first->set_facet_surface_center(f.second, p);
    }

    static const Point& get_facet_surface_center(const Facet& f)
    {
      return f.first->get_facet_surface_center(f.second);
    }

    static void reset_visited(Facet f)
    {
      f.first->reset_visited(f.second);
    }

    static bool is_facet_visited(const Facet& f)
    {
      return f.first->is_facet_visited(f.second);
    }

    // Remains unchanged
    Tr& triangulation_ref_impl()
      {
	return tr;
      }

    // Remains unchanged
    const Tr& triangulation_ref_impl() const
      {
	return tr;
      }

    // Initialization function
    void scan_triangulation_impl() {
      Point center;
      // We test only the finite Delaunay facets

#ifdef CGAL_SURFACE_MESHER_VERBOSE
      std::cout << "scanning facets..." << std::endl;
#endif
      for (Finite_facets_iterator fit = tr.finite_facets_begin(); fit !=
	     tr.finite_facets_end(); ++fit) {
	if (tr.dimension() == 3) {
          new_facet<true>(*fit);
          // see definition of
          // template <bool> new_facet()
	}
	else {
	  CGAL_assertion (tr.dimension() == 2);
	  if (is_facet_on_surface(*fit, center)) {
	    //Cell_handle c;
	    c2t3.set_in_complex(*fit);
	    set_facet_surface_center((*fit),center);
	    //c=(*fit).first;
	    //c->set_facet_on_surface((*fit).second,true);
	    //c->set_facet_surface_center((*fit).second,center);

            Quality a_r;
	    if (criteria.is_bad (*fit, a_r)) {
	      facets_to_refine.insert(*fit,a_r);
	    }
	  }
	  else
	    c2t3.remove_from_complex(*fit);
	  //	    (*fit).first->set_facet_on_surface((*fit).second,false);
	}
      }
    } // scan_triangulation_impl end

    // Tells whether there remain elements to refine
    bool no_longer_element_to_refine_impl() const {
      return facets_to_refine.empty();
    }

    // Returns the next element to refine
    Facet get_next_element_impl() {
      return facets_to_refine.front()->second;
    }

    // deletes the next element from the set of elements to refine
    // NB: it is useless here, since the update of the restricted
    // Delaunay triangulation automatically deletes the next element
    void pop_next_element_impl() {}

    // From the element to refine, gets the point to insert
    Point refinement_point_impl(const Facet& f) const
      {
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
        std::cerr << "point from Surface_mesher: ";
#endif
	CGAL_assertion (c2t3.face_status(f) == C2T3::REGULAR);
	//CGAL_assertion (f.first->is_facet_on_surface (f.second));
	return get_facet_surface_center (f);
	//return f.first->get_facet_surface_center (f.second);
      }

    // Useless here
    void before_conflicts_impl(const Facet&, const Point&
#ifdef CGAL_SURFACE_MESHER_DEBUG_BEFORE_CONFLICTS
			                                s)
    {
      std::cerr << "Refine_facets: before conflicts of " << s << " ";
#else
                                                        )
    {
#endif
    }

  // Useless here
    Mesher_level_conflict_status private_test_point_conflict_impl(const Point&,
								  Zone&)
      {
	return NO_CONFLICT;
      }

    // Useless here
    void after_no_insertion_impl(const Facet&, const Point&, const Zone&)
    {
    }

    ///////////////////////////////////////////////////////////////////////////
    // Tests if the point p encroaches one facet of the restricted Delaunay
    // triangulation.

    Mesher_level_conflict_status
    test_point_conflict_from_superior_impl(const Point& p,
					   Zone& zone)
      {
	for (typename Zone::Facets_iterator fit =
	       zone.internal_facets.begin();
	     fit != zone.internal_facets.end();
	     ++fit)
	  if( test_if_facet_is_encroached(*fit, p) )
	    return CONFLICT_BUT_ELEMENT_CAN_BE_RECONSIDERED;

	for (typename Zone::Facets_iterator fit =
	       zone.boundary_facets.begin();
	     fit != zone.boundary_facets.end();
	     ++fit)
	  if( test_if_facet_is_encroached(*fit, p) )
	    return CONFLICT_BUT_ELEMENT_CAN_BE_RECONSIDERED;

	return NO_CONFLICT;
      }

    bool test_if_facet_is_encroached(const Facet f, const Point& p)
    {
      if( tr.is_infinite(f.first) )
	return false;

      if (c2t3.face_status(f) == C2T3::REGULAR)
	{
	  const Cell_handle& c = f.first;
	  const int& index = f.second;

	  typename GT::Compute_squared_distance_3 distance =
	    tr.geom_traits().compute_squared_distance_3_object();

	  // test Delaunay surfacique
	  Point center = get_facet_surface_center(f);

          for(bool exit = false; ; exit = true)
          {
            // this for loop is a trick to pass in the following "if" once
            // with center="surface center", and once with
            // center="circumcenter"

            if( distance(center, p) <
                distance(center, c->vertex((index+1)&3)->point()) )
	    {
              Quality q;
              criteria.is_bad(f, q); // to get q (passed as reference)
              Facet other_side = mirror_facet(f);
              if(f.first < other_side.first)
                facets_to_refine.insert(f, q);
              else
                facets_to_refine.insert(other_side, q);
	      return true;
	    }

            if(exit)
              return false;

            // test Gabriel
            center = tr.geom_traits().construct_circumcenter_3_object()
              (c->vertex((index+1)&3)->point(),
               c->vertex((index+2)&3)->point(),
               c->vertex((index+3)&3)->point());
          }
	}
      return false;
    }


  /* returns the conflicts zone */
    Zone conflicts_zone_impl(const Point& p, Facet f) const {
	Zone zone;

       	// TODO may be avoid the locate here
	zone.cell =
	  tr.locate(p, zone.locate_type, zone.i, zone.j, f.first);

        tr.find_conflicts(p, zone.cell,
                          std::back_inserter(zone.boundary_facets),
                          std::back_inserter(zone.cells),
                          std::back_inserter(zone.internal_facets));
	return zone;
    }

    ///////////////////////////////////////////////////////////////////////////
    // Deletes old facets from the restricted Delaunay triangulation

    void before_insertion_impl(const Facet&, const Point& s, Zone& zone) {

      CGAL_assertion_code(Vertex_handle v);
      CGAL_assertion (!tr.is_vertex(s,v));

      if (tr.dimension() == 3) {
	// On s'occupe des facettes de la zone de conflit
	for (typename Zone::Facets_iterator fit =
	       zone.internal_facets.begin();
	     fit != zone.internal_facets.end();
	     ++fit)
	  handle_facet_inside_conflict_zone (*fit);

	for (typename Zone::Facets_iterator fit =
	       zone.boundary_facets.begin(); fit !=
	       zone.boundary_facets.end(); ++fit)
	  handle_facet_on_boundary_of_conflict_zone (*fit);
      }

      // If dim < 3, then the triangulation has only one facet and the
      // complex has no facet, generically
      else {
	CGAL_assertion (tr.dimension() == 2);
	facets_to_refine.clear();
      }

    }


    ///////////////////////////////////////////////////////////////////////////
    // Adds new facets from the restricted Delaunay triangulation

    void after_insertion_impl(const Vertex_handle& v) {