surface_mesher.h

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

#ifdef CGAL_SURFACE_MESHER_DEBUG_AFTER_INSERTION
      std::cerr << "Inserted\n";
#endif
      restore_restricted_Delaunay(v);
    }

    void restore_restricted_Delaunay(const Vertex_handle& v)
    {
      Cell_handle cellule;

      // On met a jour les flags des nouvelles facettes (celles a
      // l'interieur du trou et celles sur le bord du trou)

      std::list<Cell_handle> cellules;
      tr.incident_cells (v, std::back_inserter(cellules));

      while(!cellules.empty()) {
	cellule=cellules.front();
	cellules.pop_front();

	// Look at all four facets of the cell, starting with the
	// facet opposite to the new vertex
	int indice = cellule->index (v);
	handle_opposite_facet (Facet (cellule, indice));

	for (int i = 1; i <= 3; ++i)
	  handle_incident_facet (Facet (cellule, (indice+i)&3));
      }


    } // restore_restricted_Delaunay end




    ///////////////////////////////////////////////////////////////////////////

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




    ///////////////////////////////////////////////////////////////////////////
    // Private functions
  private:


    ///////////////////////
    // For before_insertion

    // Actions to perform on a facet inside the conflict zone
    void handle_facet_inside_conflict_zone (Facet f) {
      Facet other_side = mirror_facet(f);

      // On enleve la facette de la liste des mauvaises facettes
      if(f.first < other_side.first)
        facets_to_refine.erase(f);
      else
        facets_to_refine.erase(other_side);

      // Le compteur des visites est remis a zero
       reset_visited(f);
       reset_visited(other_side);

      // On retire la facette du complexe (car on doit etre
      // independant de l'implementation du complexe)
      c2t3.remove_from_complex (f);
    }

    // Action to perform on a facet on the boundary of the conflict zone
    void handle_facet_on_boundary_of_conflict_zone (const Facet& f) {
      // perform the same operations as for an internal facet
      handle_facet_inside_conflict_zone (f);
    }



    ///////////////////////
    // For after_insertion

    // Action to perform on a facet incident to the new vertex
    void handle_incident_facet (const Facet& f) {

      // If the facet is infinite or has been already visited,
      // then there is nothing to do as for it or its edges
      if (tr.is_infinite(f) || is_facet_visited(f))
	return;

      new_facet<false>(f);
    }

    // Action to perform on any new facet
    template <bool remove_from_complex_if_not_in_rectricted_Delaunay>
    void new_facet (const Facet& f) {
      Facet other_side = mirror_facet(f);

      // NB: set_facet_visited() is implementation dependant
      // and each side of the real facet has to be considered
      // separately
      set_facet_visited(f);
      set_facet_visited(other_side);

      // On regarde d'abord si la facette est dans le Delaunay
      // restreint
      Point center;
      if (is_facet_on_surface(f, center)) {

	// NB: set_in_complex() is implementation independant and
	// consider both sides of the facet at once
	c2t3.set_in_complex(f);

	// NB: set_facet_surface_center() is implementation dependant
	// and each side of the real facet has to be considered
	// separately
	set_facet_surface_center(f, center);
	set_facet_surface_center(other_side, center);

	// On regarde alors si la facette est bonne
        Quality a_r;
	if (criteria.is_bad (f, a_r)) {
          if(f.first < other_side.first)
            facets_to_refine.insert (f, a_r);
          else
            facets_to_refine.insert (other_side, a_r);
	}
      }
      else
        if( remove_from_complex_if_not_in_rectricted_Delaunay )
          c2t3.remove_from_complex(f);
    }

    // Action to perform on a facet opposite to the new vertex
    void handle_opposite_facet (const Facet& f) {
      // perform the same operations as for a facet incident to the
      // new vertex
      handle_incident_facet (f);
    }



    ///////////////////////
    // Predicate to test restricted Delaunay membership


    // Tests whether a given facet is restricted or not
    bool is_facet_on_surface(const Facet& f, Point& center) {
      typedef typename Surface_mesh_traits::Intersect_3 Intersect_3;
      Intersect_3 intersect = meshtraits.intersect_3_object();

      Object dual = tr.dual(f);

      //      typename GT::Segment_3 segment;
      typename GT::Segment_3 segment;
      typename GT::Ray_3 ray;
      typename GT::Line_3 line;

      Object intersection;
      // If the dual is a segment
      if (assign(segment, dual)) {
	intersection = intersect(surf, segment);
      }
      // If the dual is a ray
      else if(assign(ray, dual)) {
        // If a facet is on the convex hull, and if its finite incident
        // cell has a very bid Delaunay ball, then the dual of the facet is
        // a ray constructed with a point with very big coordinates, and a
        // vector with small coordinates. Its can happen than the
        // constructed ray is degenerate (the point(1) of the ray is
        // point(0) plus a vector whose coordinates are espilon).
        typename GT::Is_degenerate_3 is_degenerate;
        if(is_degenerate(ray)) return false;

	intersection = intersect(surf, ray);
      }
      // If the dual is a line
      else if(assign(line, dual)) {
	intersection = intersect(surf, line);
      }
      // Else there is a problem with the dual
      else {
        CGAL_assertion(false);
      }

      return assign(center,intersection);
    }

  protected:


    // Checks restricted Delaunay triangulation
    bool check_restricted_delaunay () {
      for (Finite_facets_iterator fit = tr.finite_facets_begin(); fit !=
	     tr.finite_facets_end(); ++fit) {
	Facet other_side = mirror_facet(*fit);
	CGAL_assertion (c2t3.face_status(*fit) ==
			c2t3.face_status(other_side));
	//CGAL_assertion (fit->first->is_facet_on_surface (fit->second) ==
	//		other_side.first->is_facet_on_surface
	//		(other_side.second));
	Point center;
	const bool restr = is_facet_on_surface(*fit, center);
	const bool restr_bis = is_facet_on_surface(other_side, center);
	CGAL_assertion (restr == restr_bis);
	CGAL_assertion ((c2t3.face_status(*fit)
			 == C2T3::REGULAR) == restr);
	CGAL_assertion ((c2t3.face_status(other_side)
			 == C2T3::REGULAR) == restr_bis);
	//CGAL_assertion (fit->first->is_facet_on_surface (fit->second) ==
	//		restr);
	//CGAL_assertion (other_side.first->is_facet_on_surface
	//		(other_side.second) == restr_bis);

	if ( (c2t3.face_status(*fit) == C2T3::REGULAR) !=
	    is_facet_on_surface(*fit, center)) {
	  std::cerr << "Error in restricted Delaunay triangulation: ("
		    << (c2t3.face_status(*fit) == C2T3::REGULAR)
		    << "/"
		    << is_facet_on_surface(*fit, center)
		    << ")"
		    << std::endl;
	  return false;
	}
      }
      return true;
    }

    std::string debug_info() const
    {
      std::stringstream s;
      s << facets_to_refine.size();
      return s.str();
    }

    static std::string debug_info_header()
    {
      return "number of facets";
    }
  };  // end Surface_mesher_base

    namespace details {

      template <typename Base, typename Self>
      class Mesher_level_generator {
        typedef typename Base::Complex_2_in_triangulation_3 C2T3;
        typedef typename C2T3::Triangulation Triangulation;
        typedef typename Triangulation::Facet Facet;
        typedef Triangulation_mesher_level_traits_3<Triangulation> Tr_m_l_traits_3;
      public:
        typedef Mesher_level <Triangulation,
                              Self,
                              Facet,
                              Null_mesher_level,
                              Tr_m_l_traits_3> Type;
        typedef Type type;
      }; // end class Mesher_level_generator<Base, Self>
      
    } // end namespace details

    enum Verbose_flag { VERBOSE, NOT_VERBOSE};

  template <
    typename Base,
    Verbose_flag verbose = NOT_VERBOSE
    >
  struct Surface_mesher
    : public Base,
      public details::Mesher_level_generator<Base,
                                             Surface_mesher<Base> >::Type
  {
  public:
    typedef typename Base::Complex_2_in_triangulation_3 C2T3;
    typedef typename C2T3::Triangulation Triangulation;
    typedef Triangulation Tr;
    typedef typename Base::Surface Surface;
    typedef typename Base::Criteria Criteria;
    typedef typename Base::Surface_mesh_traits Surface_mesh_traits;

    typedef Surface_mesher<Base> Self;

    typedef typename details::Mesher_level_generator<Base,
      Surface_mesher<Base> >::Type Mesher_lvl;

    using Mesher_lvl::scan_triangulation;
    using Mesher_lvl::refine;
    using Mesher_lvl::is_algorithm_done;
    using Mesher_lvl::one_step;
    using Base::check_restricted_delaunay;

    typedef C2T3 Complex_2_in_triangulation_3;

  private:
    Null_mesher_level null_mesher_level;
    Null_mesh_visitor null_visitor;
    bool initialized;

  public:
    Surface_mesher(C2T3& c2t3,
                   const Surface& surface,
                   const Surface_mesh_traits& mesh_traits,
                   const Criteria& criteria)
      : Base(c2t3, surface, mesh_traits, criteria), 
        Mesher_lvl(null_mesher_level),
        initialized(false)
    {
#ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
      std::cerr << "CONS: Surface_mesher\n";
#endif
    }

    // Initialization
    void init()
      {
	scan_triangulation();
	initialized = true;
	CGAL_assertion(check_restricted_delaunay());
      }

    void refine_mesh () {
      if(!initialized)
	init();

      if (verbose == NOT_VERBOSE)
	refine (null_visitor);
      else {
	std::cerr << "Refining...\n";
	int nbsteps = 0;
        std::cerr << "Legende of the following line: "
                  << "(number of steps," << this->debug_info_header()
                  << ")\n";
	std::cerr << "(" << nbsteps << ","
		  << this->facets_to_refine.size() << ")";
	while (!is_algorithm_done()) {
	  one_step (null_visitor);
	  std::cerr << "\r             \r"
		    << "(" << ++nbsteps << ","
		    << this->debug_info()
		    << ")";
	}
	std::cerr << "\ndone.\n";
      }

      CGAL_assertion(check_restricted_delaunay());

      initialized = false;
    }

  };  // end Surface_mesher


  }  // end namespace Surface_mesher

}  // end namespace CGAL


#endif // CGAL_SURFACE_MESHER_SURFACE_MESHER_H