1 // Copyright (c) 2013 INRIA Sophia-Antipolis (France), 2 // 2014-2015 GeometryFactory (France). 3 // All rights reserved. 4 // 5 // This file is part of CGAL (www.cgal.org). 6 // 7 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Mesh_2/include/CGAL/Mesh_2/Mesh_sizing_field.h $ 8 // $Id: Mesh_sizing_field.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 9 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 10 // 11 // 12 // Author(s) : Jane Tournois, Pierre Alliez 13 // 14 15 #ifndef CGAL_MESH_2_MESH_SIZING_FIELD_H 16 #define CGAL_MESH_2_MESH_SIZING_FIELD_H 17 18 #include <CGAL/license/Mesh_2.h> 19 20 21 #include <CGAL/Mesh_2/Sizing_field_2.h> 22 #include <CGAL/squared_distance_2.h> 23 24 namespace CGAL { 25 26 namespace Mesh_2 27 { 28 /** 29 * @class Mesh_sizing_field 30 */ 31 template <typename Tr, bool Need_vertex_update = true> 32 class Mesh_sizing_field 33 : public virtual Sizing_field_2<Tr> 34 { 35 // Types 36 typedef typename Tr::Geom_traits Gt; 37 typedef typename Tr::Point Point_2; 38 typedef typename Gt::FT FT; 39 40 typedef typename Tr::Vertex_handle Vertex_handle; 41 typedef typename Tr::Face_handle Face_handle; 42 typedef typename Tr::Edge Edge; 43 44 public: 45 // update vertices of mesh triangulation ? 46 static const bool is_vertex_update_needed = Need_vertex_update; 47 48 public: 49 /** 50 * Constructor 51 */ Mesh_sizing_field(Tr & tr)52 Mesh_sizing_field(Tr& tr) 53 : tr_(tr) 54 , last_face_() 55 { 56 init(); 57 } 58 59 /** 60 * Returns size at point \c p. 61 */ operator()62 FT operator()(const Point_2& p) const 63 { return this->operator()(p, last_face_); } 64 65 /** 66 * Returns size at point \c p, using \c v to accelerate \c p location 67 * in triangulation 68 */ operator()69 FT operator()(const Point_2& p, const Vertex_handle& v) const 70 { return this->operator()(p, v->face()); } 71 72 /** 73 * Returns size at point \c p. 74 */ operator()75 FT operator()(const Point_2& p, const Face_handle& c) const 76 { 77 const Face_handle fh = tr_.locate(p,c); 78 last_face_ = fh; 79 80 if ( !tr_.is_infinite(fh) ) 81 return interpolate_on_face_vertices(p,fh); 82 else 83 return interpolate_on_edge_vertices(p,fh); 84 } 85 86 /** 87 * Fill sizing field with actual size inside the triangulation 88 */ init()89 void init() 90 { 91 for(typename Tr::Finite_vertices_iterator 92 vit = tr_.finite_vertices_begin() ; 93 vit != tr_.finite_vertices_end() ; 94 ++vit ) 95 { 96 vit->set_sizing_info(average_incident_edge_length(vit)); 97 } 98 } 99 100 private: 101 /** 102 * Returns size at point \c p, by interpolation inside facet 103 */ interpolate_on_face_vertices(const Point_2 & p,const Face_handle & f)104 FT interpolate_on_face_vertices(const Point_2& 105 #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES 106 p 107 #endif 108 , const Face_handle& f) const 109 { 110 // Interpolate value using tet vertices values 111 const FT& sa = f->vertex(0)->sizing_info(); 112 const FT& sb = f->vertex(1)->sizing_info(); 113 const FT& sc = f->vertex(2)->sizing_info(); 114 115 #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES 116 const Point_2& a = f->vertex(0)->point(); 117 const Point_2& b = f->vertex(1)->point(); 118 const Point_2& c = f->vertex(2)->point(); 119 double abc_inv = 1. / CGAL::area(a, b, c); 120 double alpha = CGAL::area(p, b, c) * abc_inv; 121 double beta = CGAL::area(a, p, c) * abc_inv; 122 double gamma = CGAL::area(a, b, p) * abc_inv; 123 124 return alpha * sa + beta * sb + gamma * sc; 125 #else 126 return ((sa + sb + sc) / 3.); 127 #endif 128 } 129 130 /** 131 * Returns size at point \c p, by interpolation inside edge 132 * (\c e f is assumed to be an infinite face) 133 */ interpolate_on_edge_vertices(const Point_2 & p,const Face_handle & f)134 FT interpolate_on_edge_vertices(const Point_2& 135 #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES 136 p 137 #endif 138 , const Face_handle& f) const 139 { 140 CGAL_precondition(tr_.is_infinite(f)); 141 int finite_i = -1; 142 for(int i = 0; i < 3; ++i) 143 { 144 if(!tr_.is_infinite(f, i)) 145 finite_i = i; 146 } 147 CGAL_assertion(finite_i != -1); 148 149 const FT& sa = f->vertex(tr_.cw(finite_i))->sizing_info(); 150 const FT& sb = f->vertex(tr_.ccw(finite_i))->sizing_info(); 151 152 #ifdef CGAL_MESH_2_SIZING_FIELD_USE_BARYCENTRIC_COORDINATES 153 const Point_2& pa = f->vertex(tr_.cw(finite_i))->point(); 154 const Point_2& pb = f->vertex(tr_.ccw(finite_i))->point(); 155 156 FT t = CGAL::sqrt( 157 CGAL::squared_distance(pb, pb) / CGAL::squared_distance(pa, pb)); 158 CGAL_assertion(t <= 1.); 159 160 return t * sa + (1. - t) * sb; 161 #else 162 return 0.5 * (sa + sb); 163 #endif 164 } 165 average_incident_edge_length(const Vertex_handle & v)166 FT average_incident_edge_length(const Vertex_handle& v) const 167 { 168 typename Tr::Edge_circulator ec = tr_.incident_edges(v); 169 typename Tr::Edge_circulator end = ec; 170 171 FT sum_len(0.); 172 FT nb = 0.; 173 do 174 { 175 Edge e = *ec; 176 if(tr_.is_infinite(e)) 177 continue; 178 179 Face_handle f1 = e.first; 180 Face_handle f2 = e.first->neighbor(e.second); 181 if(f1->is_in_domain() || f2->is_in_domain()) 182 { 183 sum_len += length(e); 184 ++nb; 185 } 186 } 187 while(++ec != end); 188 // nb == 0 could happen if there is an isolated point. 189 if( 0 != nb ) 190 return sum_len/nb; 191 else 192 // Use outside faces to compute size of point 193 return 1.;//todo 194 } 195 length(const Edge & e)196 FT length(const Edge& e) const 197 { 198 Point_2 p1 = e.first->vertex(Tr::cw(e.second))->point(); 199 Point_2 p2 = e.first->vertex(Tr::ccw(e.second))->point(); 200 return CGAL::sqrt(CGAL::squared_distance(p1, p2)); 201 } 202 203 private: 204 /// The triangulation 205 Tr& tr_; 206 /// A face_handle that is used to accelerate location queries 207 mutable Face_handle last_face_; 208 }; 209 210 } // end namespace Mesh_2 211 212 } //namespace CGAL 213 214 #endif // CGAL_MESH_2_MESH_SIZING_FIELD_H 215