// Copyright (c) 2013 INRIA Sophia-Antipolis (France), // 2014-2015 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // You can redistribute it and/or modify it under the terms of the GNU // General Public License as published by the Free Software Foundation, // either version 3 of the License, or (at your option) any later version. // // 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$ // $Id$ // // // Author(s) : Jane Tournois, Pierre Alliez // #ifndef CGAL_MESH_2_MESH_SIZING_FIELD_H #define CGAL_MESH_2_MESH_SIZING_FIELD_H #include namespace CGAL { namespace Mesh_2 { /** * @class Mesh_sizing_field */ template class Mesh_sizing_field : public virtual Sizing_field_2 { // Types typedef typename Tr::Geom_traits Gt; typedef typename Tr::Point Point_2; typedef typename Gt::FT FT; typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Face_handle Face_handle; typedef typename Tr::Edge Edge; public: // update vertices of mesh triangulation ? static const bool is_vertex_update_needed = Need_vertex_update; public: /** * Constructor */ Mesh_sizing_field(Tr& tr) : tr_(tr) , last_face_() { init(); } /** * Returns size at point \c p. */ FT operator()(const Point_2& p) const { return this->operator()(p, last_face_); } /** * Returns size at point \c p, using \c v to accelerate \c p location * in triangulation */ FT operator()(const Point_2& p, const Vertex_handle& v) const { return this->operator()(p, v->face()); } /** * Returns size at point \c p. */ FT operator()(const Point_2& p, const Face_handle& c) const { const Face_handle fh = tr_.locate(p,c); last_face_ = fh; if ( !tr_.is_infinite(fh) ) return interpolate_on_face_vertices(p,fh); else return interpolate_on_edge_vertices(p,fh); } /** * Fill sizing field with actual size inside the triangulation */ void init() { for(typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin() ; vit != tr_.finite_vertices_end() ; ++vit ) { vit->set_sizing_info(average_incident_edge_length(vit)); } } private: /** * Returns size at point \c p, by interpolation inside facet */ FT interpolate_on_face_vertices(const Point_2& #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES p #endif , const Face_handle& f) const { // Interpolate value using tet vertices values const FT& sa = f->vertex(0)->sizing_info(); const FT& sb = f->vertex(1)->sizing_info(); const FT& sc = f->vertex(2)->sizing_info(); #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES const Point_2& a = f->vertex(0)->point(); const Point_2& b = f->vertex(1)->point(); const Point_2& c = f->vertex(2)->point(); double abc_inv = 1. / CGAL::area(a, b, c); double alpha = CGAL::area(p, b, c) * abc_inv; double beta = CGAL::area(a, p, c) * abc_inv; double gamma = CGAL::area(a, b, p) * abc_inv; return alpha * sa + beta * sb + gamma * sc; #else return ((sa + sb + sc) / 3.); #endif } /** * Returns size at point \c p, by interpolation inside edge * (\c e f is assumed to be an infinite face) */ FT interpolate_on_edge_vertices(const Point_2& #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES p #endif , const Face_handle& f) const { CGAL_precondition(tr_.is_infinite(f)); int finite_i = -1; for(int i = 0; i < 3; ++i) { if(!tr_.is_infinite(f, i)) finite_i = i; } CGAL_assertion(finite_i != -1); const FT& sa = f->vertex(tr_.cw(finite_i))->sizing_info(); const FT& sb = f->vertex(tr_.ccw(finite_i))->sizing_info(); #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES const Point_2& pa = f->vertex(tr_.cw(finite_i))->point(); const Point_2& pb = f->vertex(tr_.ccw(finite_i))->point(); FT t = CGAL::sqrt( CGAL::squared_distance(pb, pb) / CGAL::squared_distance(pa, pb)); CGAL_assertion(t <= 1.); return t * sa + (1. - t) * sb; #else return 0.5 * (sa + sb); #endif } FT average_incident_edge_length(const Vertex_handle& v) const { typename Tr::Edge_circulator ec = tr_.incident_edges(v); typename Tr::Edge_circulator end = ec; FT sum_len(0.); FT nb = 0.; do { Edge e = *ec; if(tr_.is_infinite(e)) continue; Face_handle f1 = e.first; Face_handle f2 = e.first->neighbor(e.second); if(f1->is_in_domain() || f2->is_in_domain()) { sum_len += length(e); ++nb; } } while(++ec != end); // nb == 0 could happen if there is an isolated point. if( 0 != nb ) return sum_len/nb; else // Use outside faces to compute size of point return 1.;//todo } FT length(const Edge& e) const { Point_2 p1 = e.first->vertex(Tr::cw(e.second))->point(); Point_2 p2 = e.first->vertex(Tr::ccw(e.second))->point(); return CGAL::sqrt(CGAL::squared_distance(p1, p2)); } private: /// The triangulation Tr& tr_; /// A face_handle that is used to accelerate location queries mutable Face_handle last_face_; }; } // end namespace Mesh_2 } //namespace CGAL #endif // CGAL_MESH_2_MESH_SIZING_FIELD_H