1 /* 2 * Copyright (c) 2012-2014, Bruno Levy 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions are met: 7 * 8 * * Redistributions of source code must retain the above copyright notice, 9 * this list of conditions and the following disclaimer. 10 * * Redistributions in binary form must reproduce the above copyright notice, 11 * this list of conditions and the following disclaimer in the documentation 12 * and/or other materials provided with the distribution. 13 * * Neither the name of the ALICE Project-Team nor the names of its 14 * contributors may be used to endorse or promote products derived from this 15 * software without specific prior written permission. 16 * 17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27 * POSSIBILITY OF SUCH DAMAGE. 28 * 29 * If you modify this software, you should include a notice giving the 30 * name of the person performing the modification, the date of modification, 31 * and the reason for such modification. 32 * 33 * Contact: Bruno Levy 34 * 35 * Bruno.Levy@inria.fr 36 * http://www.loria.fr/~levy 37 * 38 * ALICE Project 39 * LORIA, INRIA Lorraine, 40 * Campus Scientifique, BP 239 41 * 54506 VANDOEUVRE LES NANCY CEDEX 42 * FRANCE 43 * 44 */ 45 46 #include <geogram/mesh/mesh_degree3_vertices.h> 47 #include <geogram/mesh/mesh_geometry.h> 48 #include <geogram/mesh/index.h> 49 #include <geogram/basic/geometry_nd.h> 50 #include <geogram/basic/logger.h> 51 52 namespace { 53 54 using namespace GEO; 55 56 // TODO: include in Doxygen documentation. 57 /* 58 * There are two different numerotation for triangles: 59 * 60 * o GEO::Mesh numerotation: 61 * - Used everywhere in Vorpaline 62 * - Edge ei is right after vertex vi when turning 63 * around the facet. 64 * 65 * o Standard triangle numerotation: 66 * - Used in this file, by t_xxx() functions 67 * and Degree3Vertex class. 68 * - Edge ei is opposite to vertex vi. 69 * 70 * GEO::Mesh numerotation: 71 * v2 72 * e2 / \e1 73 * / \ 74 * v0----v1 75 * e0 76 * 77 * Standard triangle numerotation: 78 * v2 79 * e1 / \e0 80 * / \ 81 * v0----v1 82 * e2 83 */ 84 85 /** 86 * \brief Gets the vertex of a triangle by its local index, using 87 * standard triangle numerotation (edge i is opposite to vertex i). 88 * \param[in] M the mesh 89 * \param[in] t the index of the triangle 90 * \param[in] li the local index (0,1 or 2) of the vertex in \p t 91 * \return the global index of the vertex in \p M 92 * \note Everywhere else in Vorpaline, triangle edge numerotation is 93 * different from standard triangle numerotation used here. 94 */ t_vertex(const Mesh & M,index_t t,index_t li)95 inline index_t t_vertex( 96 const Mesh& M, index_t t, index_t li 97 ) { 98 geo_debug_assert(M.facets.nb_vertices(t) == 3); 99 geo_debug_assert(li < 3); 100 return M.facet_corners.vertex(M.facets.corners_begin(t) + li); 101 } 102 103 /** 104 * \brief Gets the index of an adjacent triangle, using 105 * \param[in] M the mesh 106 * \param[in] t the index of the triangle 107 * \param[in] le the local index (0,1 or 2) of the edge in \p t, using 108 * standard triangle numerotation (edge i is opposite to vertex i). 109 * \return the global index of the triangle adjacent to \p t accros 110 * edge \p le in \p M 111 * \note Everywhere else in Vorpaline, triangle edge numerotation is 112 * different from standard triangle numerotation used here. 113 */ 114 t_adjacent(const Mesh & M,index_t t,index_t le)115 inline index_t t_adjacent(const Mesh& M, index_t t, index_t le) { 116 geo_debug_assert(M.facets.nb_vertices(t) == 3); 117 geo_debug_assert(le < 3); 118 le = (le + 1) % 3; // convert corner numerotation 119 // to standard triangle numerotation 120 return M.facet_corners.adjacent_facet(M.facets.corners_begin(t) + le); 121 } 122 123 /** 124 * \brief Gets the local vertex index in a triangle by its global index, 125 * using standard triangle numerotation (edge i is opposite to vertex i). 126 * \param[in] M the mesh 127 * \param[in] t the index of the triangle 128 * \param[in] v the global index of the vertex in \p M 129 * \return the local index (0,1 or 2) of the vertex \p v in \p t 130 * \note Everywhere else in Vorpaline, triangle edge numerotation is 131 * different from standard triangle numerotation used here. 132 */ t_index(const Mesh & M,index_t t,index_t v)133 inline index_t t_index(const Mesh& M, index_t t, index_t v) { 134 geo_debug_assert(M.facets.nb_vertices(t) == 3); 135 for(index_t li = 0; li < 3; li++) { 136 if(M.facet_corners.vertex(M.facets.corners_begin(t) + li) == v) { 137 return li; 138 } 139 } 140 geo_assert_not_reached; 141 } 142 143 /** 144 * \brief Sets a triangle vertices and neighbors 145 * using standard triangle numerotation (edge i is opposite to vertex i). 146 * \param[in] M the mesh 147 * \param[in] t the index of the triangle 148 * \param[in] v1 the global index of the first vertex 149 * \param[in] v2 the global index of the second vertex 150 * \param[in] v3 the global index of the third vertex 151 * \param[in] adj1 the index of the adjacent triangle opposite to \p v1 152 * \param[in] adj2 the index of the adjacent triangle opposite to \p v2 153 * \param[in] adj3 the index of the adjacent triangle opposite to \p v3 154 * \note Everywhere else in Vorpaline, triangle edge numerotation is 155 * different from standard triangle numerotation used here. 156 */ t_set(Mesh & M,index_t t,index_t v1,index_t v2,index_t v3,index_t adj1,index_t adj2,index_t adj3)157 inline void t_set( 158 Mesh& M, 159 index_t t, 160 index_t v1, index_t v2, index_t v3, 161 index_t adj1, index_t adj2, index_t adj3 162 ) { 163 geo_debug_assert(M.facets.nb_vertices(t) == 3); 164 geo_debug_assert(adj1 != NO_FACET); 165 geo_debug_assert(adj2 != NO_FACET); 166 geo_debug_assert(adj3 != NO_FACET); 167 index_t c1 = M.facets.corners_begin(t); 168 index_t c2 = c1 + 1; 169 index_t c3 = c1 + 2; 170 M.facet_corners.set_vertex(c1, v1); 171 M.facet_corners.set_vertex(c2, v2); 172 M.facet_corners.set_vertex(c3, v3); 173 M.facet_corners.set_adjacent_facet(c1, adj3); 174 M.facet_corners.set_adjacent_facet(c2, adj1); 175 M.facet_corners.set_adjacent_facet(c3, adj2); 176 } 177 178 // TODO: ascii-art needed here !! 179 180 /** 181 * \brief Stores the three facets incident to a degree 3 vertex 182 * and the three adjacent facets. 183 * \details Implements operator< (according to distance between degree 3 184 * vertex and neighbor vertices supporting plane). 185 * \note Uses standard triangle numerotation (edge i is opposite to 186 * vertex i). Everywhere else in Vorpaline, triangle edge numerotation is 187 * different from standard triangle numerotation used here. 188 */ 189 struct Degree3Vertex { 190 191 /** 192 * \brief Constructs a new Degree3Vertex 193 * \param[in] M the mesh 194 * \param[in] v_in the index of the degree 3 vertex in \p M 195 * \param[in] t_in the index of a triangle incident to \p v_in 196 * \pre there are exactly three triangles incident to \p v_in 197 */ Degree3Vertex__anonf6c5a7d70111::Degree3Vertex198 Degree3Vertex(const Mesh& M, index_t v_in, index_t t_in) { 199 v[3] = v_in; 200 t[0] = t_in; 201 { 202 index_t i = t_index(M, t[0], v[3]); 203 index_t j = (i + 1) % 3; 204 index_t k = (j + 1) % 3; 205 206 t[1] = index_t(t_adjacent(M, t[0], j)); 207 geo_debug_assert(t[1] != index_t(-1)); 208 t[2] = index_t(t_adjacent(M, t[0], k)); 209 geo_debug_assert(t[2] != index_t(-1)); 210 v[1] = t_vertex(M, t[0], j); 211 v[2] = t_vertex(M, t[0], k); 212 adj[0] = t_adjacent(M, t[0], i); 213 } 214 215 { 216 index_t i = t_index(M, t[1], v[3]); 217 index_t j = (i + 1) % 3; 218 index_t k = (j + 1) % 3; 219 v[0] = t_vertex(M, t[1], k); 220 geo_debug_assert(t_vertex(M, t[1], j) == v[2]); 221 geo_debug_assert(t_adjacent(M, t[1], j) == t[2]); 222 geo_debug_assert(t_adjacent(M, t[1], k) == t[0]); 223 adj[1] = t_adjacent(M, t[1], i); 224 } 225 226 { 227 index_t i = t_index(M, t[2], v[3]); 228 #ifdef GEO_DEBUG 229 index_t j = (i + 1) % 3; 230 index_t k = (j + 1) % 3; 231 geo_debug_assert(t_vertex(M, t[2], j) == v[0]); 232 geo_debug_assert(t_vertex(M, t[2], k) == v[1]); 233 geo_debug_assert(t_adjacent(M, t[2], j) == t[0]); 234 geo_debug_assert(t_adjacent(M, t[2], k) == t[1]); 235 #endif 236 adj[2] = t_adjacent(M, t[2], i); 237 } 238 239 const vec3& p0 = Geom::mesh_vertex(M, v[0]); 240 const vec3& p1 = Geom::mesh_vertex(M, v[1]); 241 const vec3& p2 = Geom::mesh_vertex(M, v[2]); 242 const vec3& p3 = Geom::mesh_vertex(M, v[3]); 243 244 dist = ::sqrt( 245 Geom::point_triangle_squared_distance(p3, p0, p1, p2) 246 ); 247 } 248 249 /** 250 * \brief Compares two Degree3Vertex by the distance to the 251 * supporting planes of the neighbors. 252 * \param[in] rhs the comparand 253 * \return true if this Degree3Vertex can be removed before 254 * \p rhs, false otherwise 255 */ operator <__anonf6c5a7d70111::Degree3Vertex256 bool operator< (const Degree3Vertex& rhs) const { 257 return dist < rhs.dist; 258 } 259 260 double dist; 261 index_t v[4]; 262 index_t t[3]; 263 index_t adj[3]; 264 }; 265 } 266 267 /****************************************************************************/ 268 269 namespace GEO { 270 remove_degree3_vertices(Mesh & M,double max_dist)271 index_t remove_degree3_vertices(Mesh& M, double max_dist) { 272 273 if(max_dist == 0.0) { 274 return 0; 275 } 276 277 // Step 1: detect degree3 vertices 278 279 vector<signed_index_t> vertex_degree(M.vertices.nb(), 0); 280 // or -1 if v is on border or if v has an incident facet that 281 // is not a triangle. 282 283 for(index_t f: M.facets) { 284 bool f_is_triangle = (M.facets.nb_vertices(f) == 3); 285 for(index_t c: M.facets.corners(f)) { 286 index_t v = M.facet_corners.vertex(c); 287 if( 288 !f_is_triangle || 289 (M.facet_corners.adjacent_facet(c) == NO_FACET) 290 ) { 291 vertex_degree[v] = -1; 292 } else { 293 if(vertex_degree[v] != -1) { 294 vertex_degree[v]++; 295 } 296 } 297 } 298 } 299 300 // Step 2: count degree3 vertices 301 302 index_t nb_degree3_vertices = 0; 303 for(index_t v: M.vertices) { 304 if(vertex_degree[v] == 3) { 305 nb_degree3_vertices++; 306 } 307 } 308 309 if(nb_degree3_vertices == 0) { 310 Logger::out("Degree3") 311 << "Does not have any degree 3 vertex (good)" << std::endl; 312 return 0; 313 } 314 315 // Step 3: v2f[v] is one of the facets adjacent to v 316 317 vector<index_t> v2f(M.vertices.nb()); 318 for(index_t f: M.facets) { 319 for(index_t c: M.facets.corners(f)) { 320 index_t v = M.facet_corners.vertex(c); 321 v2f[v] = f; 322 } 323 } 324 325 // Step 4: compute and sort degree3 vertices configurations 326 327 vector<Degree3Vertex> degree3vertices; 328 degree3vertices.reserve(nb_degree3_vertices); 329 for(index_t v: M.vertices) { 330 if(vertex_degree[v] == 3) { 331 Degree3Vertex V(M, v, v2f[v]); 332 if(V.dist < max_dist) { 333 degree3vertices.push_back(V); 334 } 335 } 336 } 337 338 Logger::out("Degree3") 339 << "Removing " << degree3vertices.size() 340 << "/" << nb_degree3_vertices 341 << " degree 3 vertices (within max_deg3_dist)" 342 << std::endl; 343 344 std::sort(degree3vertices.begin(), degree3vertices.end()); 345 346 // Step 5: remove degree3 vertices in order 347 348 enum FacetStatus { 349 F_UNTOUCHED = 0, 350 F_TO_REMOVE = 1, 351 F_MODIFIED = 2 352 }; 353 354 vector<index_t> facet_status(M.facets.nb(), F_UNTOUCHED); 355 index_t nb_removed = 0; 356 357 for(index_t i = 0; i < degree3vertices.size(); i++) { 358 const Degree3Vertex& V = degree3vertices[i]; 359 index_t t1 = V.t[0]; 360 index_t t2 = V.t[1]; 361 index_t t3 = V.t[2]; 362 363 // Skip vertex if one of its incident facet was already modified 364 // when removing a neighboring degree 3 vertex. 365 if( 366 facet_status[t1] != F_UNTOUCHED || 367 facet_status[t2] != F_UNTOUCHED || 368 facet_status[t3] != F_UNTOUCHED 369 ) { 370 continue; 371 } 372 373 // Skip vertex if one of its neighboring facet 374 // was already modified... 375 index_t adj1 = V.adj[0]; 376 index_t adj2 = V.adj[1]; 377 index_t adj3 = V.adj[2]; 378 if(adj1 != NO_FACET && facet_status[adj1] != F_UNTOUCHED) { 379 continue; 380 } 381 if(adj2 != NO_FACET && facet_status[adj2] != F_UNTOUCHED) { 382 continue; 383 } 384 if(adj3 != NO_FACET && facet_status[adj3] != F_UNTOUCHED) { 385 continue; 386 } 387 388 facet_status[t2] = F_TO_REMOVE; // t2 will be deleted. 389 facet_status[t3] = F_TO_REMOVE; // t3 will be deleted. 390 facet_status[t1] = F_MODIFIED; // t1 is recycled and 391 // used by the new facet. 392 t_set( 393 M, t1, V.v[0], V.v[1], V.v[2], 394 index_t(V.adj[0]), index_t(V.adj[1]), index_t(V.adj[2]) 395 ); 396 397 // connect the neighbors of t2 and t3 to t1 398 for(index_t j = 1; j <= 2; j++) { 399 if(V.adj[j] != NO_FACET) { 400 index_t f = index_t(V.adj[j]); 401 for(index_t c: M.facets.corners(f)) { 402 if(M.facet_corners.adjacent_facet(c) == V.t[j]) { 403 M.facet_corners.set_adjacent_facet(c, V.t[0]); 404 } 405 } 406 } 407 } 408 nb_removed++; 409 } 410 411 // Step 6: delete dangling triangles 412 413 // Do not delete the facets that were recycled. 414 for(index_t i = 0; i < facet_status.size(); i++) { 415 if(facet_status[i] == F_MODIFIED) { 416 facet_status[i] = 0; 417 } 418 } 419 420 // Note: remove_facets() calls remove_isolated_vertices() 421 M.facets.delete_elements(facet_status); 422 423 return nb_removed; 424 } 425 } 426 427