| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2022 Inria | ||
| 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 | * Contact: Bruno Levy | ||
| 30 | * | ||
| 31 | * https://www.inria.fr/fr/bruno-levy | ||
| 32 | * | ||
| 33 | * Inria, | ||
| 34 | * Domaine de Voluceau, | ||
| 35 | * 78150 Le Chesnay - Rocquencourt | ||
| 36 | * FRANCE | ||
| 37 | * | ||
| 38 | */ | ||
| 39 | |||
| 40 | #include <geogram/voronoi/RVD_callback.h> | ||
| 41 | #include <geogram/voronoi/RVD_mesh_builder.h> | ||
| 42 | #include <geogram/voronoi/generic_RVD_cell.h> | ||
| 43 | #include <geogram/mesh/mesh_geometry.h> | ||
| 44 | #include <geogram/mesh/mesh_io.h> | ||
| 45 | #include <geogram/basic/argused.h> | ||
| 46 | |||
| 47 | namespace { | ||
| 48 | using namespace GEO; | ||
| 49 | |||
| 50 | /** | ||
| 51 | * \brief Gets the maximum region index in the mesh, not counting | ||
| 52 | * index_t(-1) triangles. | ||
| 53 | * \return The maximum region index, or index_t(-1) if all region | ||
| 54 | * indices are index_t(-1) | ||
| 55 | */ | ||
| 56 | ✗ | index_t max_region(Mesh& mesh, Attribute<index_t>& facet_region) { | |
| 57 | ✗ | index_t result = index_t(-1); | |
| 58 | ✗ | for(index_t f=0; f<mesh.facets.nb(); ++f) { | |
| 59 | ✗ | if(facet_region[f] != index_t(-1)) { | |
| 60 | ✗ | if(result == index_t(-1)) { | |
| 61 | ✗ | result = facet_region[f]; | |
| 62 | } else { | ||
| 63 | ✗ | result = std::max(result, facet_region[f]); | |
| 64 | } | ||
| 65 | } | ||
| 66 | } | ||
| 67 | ✗ | return result; | |
| 68 | } | ||
| 69 | |||
| 70 | /** | ||
| 71 | * \brief Splits the regions along hard edges. | ||
| 72 | * \param[in,out] mesh a reference to the surface mesh. | ||
| 73 | * \param[in,out] facet_region the attribute that defines the regions. | ||
| 74 | * The regions with facet_region == index_t(-1) will be split. | ||
| 75 | * \param[in] threshold (in degrees). Edges with two adjacent facets | ||
| 76 | * with normals that form an angle larger than \p threshold degrees are | ||
| 77 | * considered as hard edges. | ||
| 78 | */ | ||
| 79 | ✗ | void split_regions_along_hard_edges( | |
| 80 | Mesh& mesh, Attribute<index_t>& facet_region, double threshold | ||
| 81 | ) { | ||
| 82 | |||
| 83 | ✗ | threshold *= (M_PI / 180.0); | |
| 84 | |||
| 85 | ✗ | index_t R = max_region(mesh, facet_region); | |
| 86 | ✗ | if(R == index_t(-1)) { | |
| 87 | ✗ | R = 0; | |
| 88 | } else { | ||
| 89 | ✗ | ++R; | |
| 90 | } | ||
| 91 | ✗ | vector<bool> is_crease(mesh.facet_corners.nb(), false); | |
| 92 | ✗ | for(index_t f=0; f<mesh.facets.nb(); ++f) { | |
| 93 | ✗ | vec3 n = Geom::mesh_facet_normal(mesh,f); | |
| 94 | ✗ | for( | |
| 95 | ✗ | index_t c=mesh.facets.corners_begin(f); | |
| 96 | ✗ | c<mesh.facets.corners_end(f); ++c | |
| 97 | ) { | ||
| 98 | ✗ | index_t f2 = mesh.facet_corners.adjacent_facet(c); | |
| 99 | ✗ | if(f2 != index_t(-1)) { | |
| 100 | ✗ | vec3 n2 = Geom::mesh_facet_normal(mesh,f2); | |
| 101 | ✗ | double alpha = Geom::angle(n,n2); | |
| 102 | ✗ | if(alpha > threshold) { | |
| 103 | ✗ | is_crease[c] = true; | |
| 104 | } | ||
| 105 | } | ||
| 106 | } | ||
| 107 | } | ||
| 108 | ✗ | for(index_t f=0; f<mesh.facets.nb(); ++f) { | |
| 109 | ✗ | if(facet_region[f] == index_t(-1)) { | |
| 110 | ✗ | std::stack<index_t> S; | |
| 111 | ✗ | facet_region[f] = R; | |
| 112 | ✗ | S.push(f); | |
| 113 | ✗ | while(!S.empty()) { | |
| 114 | ✗ | index_t f2 = S.top(); | |
| 115 | ✗ | S.pop(); | |
| 116 | ✗ | for( | |
| 117 | ✗ | index_t c=mesh.facets.corners_begin(f2); | |
| 118 | ✗ | c<mesh.facets.corners_end(f2); ++c | |
| 119 | ) { | ||
| 120 | ✗ | if(is_crease[c]) { | |
| 121 | ✗ | continue; | |
| 122 | } | ||
| 123 | ✗ | index_t f3 = mesh.facet_corners.adjacent_facet(c); | |
| 124 | ✗ | if( | |
| 125 | ✗ | f3 != index_t(-1) && | |
| 126 | ✗ | facet_region[f3] == index_t(-1) | |
| 127 | ) { | ||
| 128 | ✗ | facet_region[f3] = R; | |
| 129 | ✗ | S.push(f3); | |
| 130 | } | ||
| 131 | } | ||
| 132 | } | ||
| 133 | ✗ | ++R; | |
| 134 | ✗ | } | |
| 135 | } | ||
| 136 | ✗ | } | |
| 137 | |||
| 138 | /** | ||
| 139 | * \brief Simplifies the facets of a surface mesh based on an attribute. | ||
| 140 | * \details Groups of connected facets with the same attribute value are | ||
| 141 | * replaced with a single facet. | ||
| 142 | * \param[in,out] mesh a reference to the surface mesh to be simplified | ||
| 143 | * \param[in] facet_region a reference to the attribute with the facet | ||
| 144 | * region | ||
| 145 | * \param[in] angle_threshold (in degrees). In the outer region | ||
| 146 | * (i.e. facet_region == index_t(-1)), an edge shared by two adjacent | ||
| 147 | * facets is suppressed if the angle between the | ||
| 148 | * two facet normals is smaller than \p angle_threshold. | ||
| 149 | */ | ||
| 150 | ✗ | bool simplify( | |
| 151 | Mesh& mesh, | ||
| 152 | Attribute<index_t>& facet_region, | ||
| 153 | double angle_threshold | ||
| 154 | ) { | ||
| 155 | |||
| 156 | ✗ | bool keep_outer_region = (angle_threshold == 0.0); | |
| 157 | |||
| 158 | ✗ | index_t max_r = 0; | |
| 159 | |||
| 160 | ✗ | if(!keep_outer_region) { | |
| 161 | ✗ | max_r = max_region(mesh, facet_region); | |
| 162 | ✗ | split_regions_along_hard_edges(mesh, facet_region, angle_threshold); | |
| 163 | } | ||
| 164 | |||
| 165 | ✗ | vector<bool> is_corner(mesh.vertices.nb(),false); | |
| 166 | |||
| 167 | { | ||
| 168 | ✗ | vector<index_t> rgn1(mesh.vertices.nb(), index_t(-2)); | |
| 169 | ✗ | vector<index_t> rgn2(mesh.vertices.nb(), index_t(-2)); | |
| 170 | |||
| 171 | // Keep all the vertices adjacent to 3 regions or more | ||
| 172 | ✗ | for(index_t f=0; f<mesh.facets.nb(); ++f) { | |
| 173 | ✗ | index_t r = facet_region[f]; | |
| 174 | ✗ | for(index_t lv=0; lv<mesh.facets.nb_vertices(f); ++lv) { | |
| 175 | ✗ | index_t v = mesh.facets.vertex(f,lv); | |
| 176 | ✗ | if(rgn1[v] == r || rgn2[v] == r) { | |
| 177 | ✗ | continue; | |
| 178 | } | ||
| 179 | ✗ | if(rgn1[v] == index_t(-2)) { | |
| 180 | ✗ | rgn1[v] = r; | |
| 181 | ✗ | } else if(rgn2[v] == index_t(-2)) { | |
| 182 | ✗ | rgn2[v] = r; | |
| 183 | } else { | ||
| 184 | ✗ | is_corner[v] = true; | |
| 185 | } | ||
| 186 | } | ||
| 187 | } | ||
| 188 | |||
| 189 | // Keep also the vertices adjacent to region -1 and | ||
| 190 | // to another region | ||
| 191 | |||
| 192 | ✗ | if(keep_outer_region) { | |
| 193 | ✗ | for(index_t v=0; v<mesh.vertices.nb(); ++v) { | |
| 194 | ✗ | if( | |
| 195 | ✗ | (rgn1[v] == index_t(-1) && rgn2[v] != index_t(-1)) || | |
| 196 | ✗ | (rgn2[v] == index_t(-1) && rgn1[v] != index_t(-1)) | |
| 197 | ) { | ||
| 198 | ✗ | is_corner[v] = true; | |
| 199 | } | ||
| 200 | } | ||
| 201 | } | ||
| 202 | ✗ | } | |
| 203 | |||
| 204 | // -1: not visited, 0: keep, 1: delete | ||
| 205 | ✗ | vector<index_t> facet_status(mesh.facets.nb(),index_t(-1)); | |
| 206 | |||
| 207 | // Needs to be backed-up, we are modifying the mesh !! | ||
| 208 | ✗ | index_t nf = mesh.facets.nb(); | |
| 209 | |||
| 210 | ✗ | for(index_t f=0; f<nf; ++f) { | |
| 211 | ✗ | if(facet_status[f] == index_t(-1)) { | |
| 212 | ✗ | index_t r = facet_region[f]; | |
| 213 | ✗ | if(keep_outer_region && (r == index_t(-1))) { | |
| 214 | ✗ | facet_status[f] = 0; | |
| 215 | ✗ | continue; | |
| 216 | } else { | ||
| 217 | ✗ | std::stack<index_t> S; | |
| 218 | ✗ | std::map<index_t, index_t> border_next; | |
| 219 | ✗ | S.push(f); | |
| 220 | ✗ | facet_status[f] = 1; | |
| 221 | ✗ | while(!S.empty()) { | |
| 222 | ✗ | index_t f2 = S.top(); | |
| 223 | ✗ | S.pop(); | |
| 224 | ✗ | for(index_t c1=mesh.facets.corners_begin(f2); | |
| 225 | ✗ | c1<mesh.facets.corners_end(f2); ++c1) { | |
| 226 | ✗ | index_t f3 = mesh.facet_corners.adjacent_facet(c1); | |
| 227 | ✗ | if(f3 == index_t(-1) || facet_region[f3] != r) { | |
| 228 | index_t c2 = | ||
| 229 | ✗ | mesh.facets.next_corner_around_facet(f2,c1); | |
| 230 | ✗ | index_t v1 = mesh.facet_corners.vertex(c1); | |
| 231 | ✗ | index_t v2 = mesh.facet_corners.vertex(c2); | |
| 232 | ✗ | if(border_next.find(v1) != border_next.end()) { | |
| 233 | ✗ | Logger::warn("Simplify") | |
| 234 | ✗ | << "Region has non-manifold border" | |
| 235 | ✗ | << std::endl; | |
| 236 | // Yes, goto!, why not, what's wrong with goto ? | ||
| 237 | // Why would it be ok to throw() and not goto ? | ||
| 238 | ✗ | goto rollback; | |
| 239 | } | ||
| 240 | ✗ | border_next[v1] = v2; | |
| 241 | } else { | ||
| 242 | ✗ | if(facet_status[f3] == index_t(-1)) { | |
| 243 | ✗ | facet_status[f3] = 1; | |
| 244 | ✗ | S.push(f3); | |
| 245 | } | ||
| 246 | } | ||
| 247 | } | ||
| 248 | } | ||
| 249 | ✗ | index_t nb_border_visited = 0; | |
| 250 | ✗ | vector<index_t> new_facet; | |
| 251 | ✗ | index_t v = border_next.begin()->first; | |
| 252 | do { | ||
| 253 | ✗ | if(is_corner[v]) { | |
| 254 | ✗ | new_facet.push_back(v); | |
| 255 | } | ||
| 256 | ✗ | ++nb_border_visited; | |
| 257 | ✗ | v = border_next[v]; | |
| 258 | ✗ | if(nb_border_visited > mesh.vertices.nb()) { | |
| 259 | ✗ | Logger::warn("Simplify") | |
| 260 | ✗ | << "Region has singular border topology" | |
| 261 | ✗ | << std::endl; | |
| 262 | ✗ | goto rollback; | |
| 263 | } | ||
| 264 | ✗ | } while(v != border_next.begin()->first); | |
| 265 | |||
| 266 | ✗ | if(nb_border_visited != border_next.size()) { | |
| 267 | ✗ | Logger::warn("Simplify") | |
| 268 | ✗ | << "Region has multiple borders" | |
| 269 | ✗ | << std::endl; | |
| 270 | ✗ | goto rollback; | |
| 271 | } | ||
| 272 | |||
| 273 | ✗ | if(new_facet.size() < 3) { | |
| 274 | ✗ | Logger::warn("Simplify") | |
| 275 | ✗ | << "Region has border with less than 3 corners" | |
| 276 | ✗ | << std::endl; | |
| 277 | ✗ | goto rollback; | |
| 278 | } | ||
| 279 | |||
| 280 | ✗ | index_t new_f = mesh.facets.nb(); | |
| 281 | ✗ | mesh.facets.create_polygon(new_facet.size()); | |
| 282 | ✗ | for(index_t i=0; i<new_facet.size(); ++i) { | |
| 283 | ✗ | mesh.facets.set_vertex(new_f,i,new_facet[i]); | |
| 284 | } | ||
| 285 | ✗ | facet_region[new_f] = r; | |
| 286 | ✗ | } | |
| 287 | } | ||
| 288 | } | ||
| 289 | |||
| 290 | ✗ | facet_status.resize(mesh.facets.nb(), 0); | |
| 291 | ✗ | mesh.facets.delete_elements(facet_status); | |
| 292 | ✗ | if(!keep_outer_region) { | |
| 293 | ✗ | for(index_t f=0; f<mesh.facets.nb(); ++f) { | |
| 294 | ✗ | if(facet_region[f] > max_r) { | |
| 295 | ✗ | facet_region[f] = index_t(-1); | |
| 296 | } | ||
| 297 | } | ||
| 298 | } | ||
| 299 | ✗ | return true; | |
| 300 | |||
| 301 | ✗ | rollback: | |
| 302 | ✗ | Logger::out("Simplify") << "...Rolling back." << std::endl; | |
| 303 | // delete all facets... | ||
| 304 | ✗ | facet_status.resize(mesh.facets.nb(), 1); | |
| 305 | // ... except the initial ones ! | ||
| 306 | ✗ | for(index_t f=0; f<nf; ++f) { | |
| 307 | ✗ | facet_status[f] = 0; | |
| 308 | } | ||
| 309 | ✗ | mesh.facets.delete_elements(facet_status); | |
| 310 | ✗ | if(!keep_outer_region) { | |
| 311 | ✗ | for(index_t f=0; f<mesh.facets.nb(); ++f) { | |
| 312 | ✗ | if(facet_region[f] > max_r) { | |
| 313 | ✗ | facet_region[f] = index_t(-1); | |
| 314 | } | ||
| 315 | } | ||
| 316 | } | ||
| 317 | ✗ | return false; | |
| 318 | ✗ | } | |
| 319 | |||
| 320 | /** | ||
| 321 | * \brief Gets a 2d polygon that represents a mesh facet. | ||
| 322 | * \param[in] mesh a const reference to the mesh | ||
| 323 | * \param[in] f the facet | ||
| 324 | * \param[in] N the normal vector to the facet | ||
| 325 | * \param[out] P the vertices of the polygon | ||
| 326 | * \param[out] P_ind the global indices of the vertices in \p mesh | ||
| 327 | */ | ||
| 328 | ✗ | void get_mesh_polygon2d( | |
| 329 | const Mesh& mesh, | ||
| 330 | index_t f, | ||
| 331 | const vec3& N, | ||
| 332 | vector<vec2>& P, | ||
| 333 | vector<index_t>& P_ind | ||
| 334 | ) { | ||
| 335 | ✗ | P.resize(0); | |
| 336 | ✗ | P_ind.resize(0); | |
| 337 | ✗ | vec3 Z = normalize(N); | |
| 338 | ✗ | vec3 X = Geom::perpendicular(Z); | |
| 339 | ✗ | vec3 Y = cross(Z,X); | |
| 340 | ✗ | vec3 C = Geom::mesh_facet_center(mesh,f); | |
| 341 | ✗ | index_t n = mesh.facets.nb_vertices(f); | |
| 342 | ✗ | FOR(lv,n) { | |
| 343 | ✗ | index_t v = mesh.facets.vertex(f,lv); | |
| 344 | ✗ | vec3 W = mesh.vertices.point(v)-C; | |
| 345 | ✗ | P_ind.push_back(v); | |
| 346 | ✗ | P.push_back(vec2(dot(W,X), dot(W,Y))); | |
| 347 | } | ||
| 348 | // TODO: normalize vertices order so that two | ||
| 349 | // opposite facets will have the same tessellation. | ||
| 350 | ✗ | } | |
| 351 | |||
| 352 | |||
| 353 | /** | ||
| 354 | * \brief Evaluates the score of a triangle in a closed polygon. | ||
| 355 | * \param[in] pts the closed polygon | ||
| 356 | * \param[in] i , j , k the three vertices of the triangle | ||
| 357 | * \retval 1024 if a concave angle was encountered or if the proposed | ||
| 358 | * triangle contains one of the points. | ||
| 359 | * \retval the maximum angle of the proposed triangle otherwise. | ||
| 360 | */ | ||
| 361 | ✗ | double triangle_cost( | |
| 362 | const vector<vec2>& pts, index_t i, index_t j, index_t k | ||
| 363 | ) { | ||
| 364 | ✗ | vec2 C[3] = { pts[i], pts[j], pts[k] }; | |
| 365 | ✗ | double m = 0; | |
| 366 | ✗ | FOR(v, 3) { | |
| 367 | // note that angle is not the angle inside the triangle, | ||
| 368 | // but its complement | ||
| 369 | // angle variable has the "direction" information, thus it | ||
| 370 | // is negative for concave angles (right turn) and positive | ||
| 371 | // for convex angles (left turn) | ||
| 372 | ✗ | double angle = atan2( | |
| 373 | det( | ||
| 374 | ✗ | C[(v + 1) % 3] - C[(v + 0) % 3], | |
| 375 | ✗ | C[(v + 2) % 3] - C[(v + 1) % 3] | |
| 376 | ), | ||
| 377 | dot( | ||
| 378 | ✗ | C[(v + 1) % 3] - C[(v + 0) % 3], | |
| 379 | ✗ | C[(v + 2) % 3] - C[(v + 1) % 3] | |
| 380 | ) | ||
| 381 | ); | ||
| 382 | ✗ | if (angle <= 0) return 1024.; | |
| 383 | ✗ | m = std::max(m, M_PI - angle); | |
| 384 | } | ||
| 385 | |||
| 386 | ✗ | FOR(other, pts.size()) { | |
| 387 | // TODO: check also whether triangle is inversed ? | ||
| 388 | // To be checked: I think it is already done in | ||
| 389 | // angle computations above. | ||
| 390 | ✗ | if (other == i || other == j || other == k) { | |
| 391 | ✗ | continue; | |
| 392 | } | ||
| 393 | ✗ | const vec2& P = pts[other]; | |
| 394 | ✗ | bool inside = true; | |
| 395 | ✗ | FOR(l, 3) { | |
| 396 | ✗ | inside = inside && (det(C[(l + 1) % 3] - C[l], P - C[l]) > 0); | |
| 397 | } | ||
| 398 | ✗ | if (inside) { | |
| 399 | ✗ | return 1024.0; | |
| 400 | } | ||
| 401 | } | ||
| 402 | ✗ | return m; | |
| 403 | } | ||
| 404 | |||
| 405 | /** | ||
| 406 | * \brief Triangulates a (possibly non-convex) polygon. | ||
| 407 | * \note The algorithm is in O(n^4) (bad but good enough for now). | ||
| 408 | * \param[in] pts the polygon | ||
| 409 | * \param[out] triangles the indices of the triangles vertices | ||
| 410 | * \retval true on success | ||
| 411 | * \retval false otherwise | ||
| 412 | */ | ||
| 413 | ✗ | bool triangulate_polygon( | |
| 414 | const vector<vec2>& pts, vector<index_t>& triangles | ||
| 415 | ) { | ||
| 416 | ✗ | triangles.resize(0); | |
| 417 | ✗ | index_t n = pts.size(); | |
| 418 | ✗ | geo_assert(n >= 3); | |
| 419 | |||
| 420 | ✗ | if (n == 3) { | |
| 421 | ✗ | FOR(v, 3) { | |
| 422 | ✗ | triangles.push_back(v); | |
| 423 | } | ||
| 424 | ✗ | return true; | |
| 425 | } | ||
| 426 | |||
| 427 | // we store in this table results of subproblems | ||
| 428 | // table[i*n + j] stores the triangulation cost for points from i to j | ||
| 429 | // the entry table[0*n + n-1] has the final result. | ||
| 430 | ✗ | vector<double> table(n*n, 0.); | |
| 431 | |||
| 432 | // this table stores triangle indices: | ||
| 433 | // for each subproblem (i,j) we have table[i*n + j]==k, | ||
| 434 | // i.e. the triangle is (i,k,j) | ||
| 435 | ✗ | vector<index_t> tri(n*n, index_t(-1)); | |
| 436 | |||
| 437 | // note that the table is filled in diagonals; | ||
| 438 | // elements below main diagonal are not used at all | ||
| 439 | ✗ | for (index_t pbsize = 2; pbsize < n; pbsize++) { | |
| 440 | ✗ | for (index_t i = 0, j = pbsize; j < n; i++, j++) { | |
| 441 | // recall that we are testing triangle (i,k,j) | ||
| 442 | // which splits the problem (i,j) into | ||
| 443 | // two smaller subproblems (i,k) and (k,j) | ||
| 444 | |||
| 445 | ✗ | double minv = 1e20; | |
| 446 | |||
| 447 | ✗ | index_t mink = index_t(-1); | |
| 448 | |||
| 449 | ✗ | for (index_t k = i + 1; k < j; k++) { | |
| 450 | |||
| 451 | double val = | ||
| 452 | ✗ | table[i*n + k] + table[k*n + j] + | |
| 453 | ✗ | triangle_cost(pts, i, k, j); | |
| 454 | |||
| 455 | ✗ | if (minv <= val) { | |
| 456 | ✗ | continue; | |
| 457 | } | ||
| 458 | ✗ | minv = val; | |
| 459 | ✗ | mink = k; | |
| 460 | } | ||
| 461 | ✗ | geo_assert(mink!=index_t(-1)); | |
| 462 | ✗ | table[i*n + j] = minv; | |
| 463 | ✗ | tri[i*n + j] = mink; | |
| 464 | } | ||
| 465 | } | ||
| 466 | |||
| 467 | ✗ | vector<index_t> Q(1, n - 1); | |
| 468 | ✗ | FOR(t, Q.size()) { | |
| 469 | ✗ | index_t idx = Q[t]; | |
| 470 | |||
| 471 | ✗ | index_t i = idx / n; | |
| 472 | ✗ | index_t k = tri[idx]; | |
| 473 | ✗ | index_t j = idx % n; | |
| 474 | |||
| 475 | ✗ | geo_assert(i!=index_t(-1) && k != index_t(-1) && j!=index_t(-1)); | |
| 476 | |||
| 477 | ✗ | triangles.push_back(i); | |
| 478 | ✗ | triangles.push_back(k); | |
| 479 | ✗ | triangles.push_back(j); | |
| 480 | |||
| 481 | ✗ | if (k + 2 <= j) { | |
| 482 | ✗ | Q.push_back(k*n + j); | |
| 483 | } | ||
| 484 | ✗ | if (i + 2 <= k) { | |
| 485 | ✗ | Q.push_back(i*n + k); | |
| 486 | } | ||
| 487 | } | ||
| 488 | ✗ | return table[n-1] < 1024.; | |
| 489 | ✗ | } | |
| 490 | |||
| 491 | /** | ||
| 492 | * \brief Tests whether a 2d polygon is convex. | ||
| 493 | * \param[in] P a const reference to the polygon. | ||
| 494 | * \retval true if the polygon \p P is convex. | ||
| 495 | * \retval false otherwise. | ||
| 496 | */ | ||
| 497 | ✗ | bool polygon_is_convex(const vector<vec2>& P) { | |
| 498 | ✗ | Sign s = ZERO; | |
| 499 | ✗ | FOR(i, P.size()) { | |
| 500 | ✗ | index_t j = (i+1)%P.size(); | |
| 501 | ✗ | index_t k = (j+1)%P.size(); | |
| 502 | ✗ | Sign cur_s = PCK::orient_2d(P[i], P[j], P[k]); | |
| 503 | ✗ | if(int(cur_s) * int(s) == -1) { | |
| 504 | ✗ | return false; | |
| 505 | } | ||
| 506 | ✗ | if(s == ZERO) { | |
| 507 | ✗ | s = cur_s; | |
| 508 | } | ||
| 509 | } | ||
| 510 | ✗ | return true; | |
| 511 | } | ||
| 512 | |||
| 513 | /** | ||
| 514 | * \brief Tesselates the non-convex facets of a mesh. | ||
| 515 | * \param[in,out] mesh a pointer to the mesh. | ||
| 516 | */ | ||
| 517 | ✗ | void tessellate_non_convex_facets( | |
| 518 | Mesh* mesh | ||
| 519 | ) { | ||
| 520 | // TODO: use facet_seed_ attribute and replace normal vector | ||
| 521 | // with (seed-facet seed) vector. | ||
| 522 | ✗ | vector<index_t> to_delete; | |
| 523 | ✗ | vector<vec2> P; | |
| 524 | ✗ | vector<index_t> P_ind; | |
| 525 | ✗ | vector<index_t> P_tri; | |
| 526 | ✗ | index_t nf = mesh->facets.nb(); | |
| 527 | ✗ | FOR(f, nf) { | |
| 528 | ✗ | vec3 N = Geom::mesh_facet_normal(*mesh, f); | |
| 529 | ✗ | get_mesh_polygon2d(*mesh, f, N, P, P_ind); | |
| 530 | ✗ | if(!polygon_is_convex(P)) { | |
| 531 | ✗ | if(triangulate_polygon(P, P_tri)) { | |
| 532 | ✗ | to_delete.resize(mesh->facets.nb(),0); | |
| 533 | ✗ | to_delete[f] = 1; | |
| 534 | ✗ | FOR(t, P_tri.size()/3) { | |
| 535 | ✗ | index_t newf = mesh->facets.create_triangle( | |
| 536 | ✗ | P_ind[P_tri[3*t ]], | |
| 537 | ✗ | P_ind[P_tri[3*t+1]], | |
| 538 | ✗ | P_ind[P_tri[3*t+2]] | |
| 539 | ); | ||
| 540 | ✗ | mesh->facets.attributes().copy_item(newf,f); | |
| 541 | } | ||
| 542 | } else { | ||
| 543 | ✗ | Logger::warn("RVD") | |
| 544 | ✗ | << "Could not triangulate non-convex facet" | |
| 545 | ✗ | << std::endl; | |
| 546 | } | ||
| 547 | } | ||
| 548 | } | ||
| 549 | ✗ | if(to_delete.size() != 0) { | |
| 550 | ✗ | to_delete.resize(mesh->facets.nb(), 0); | |
| 551 | ✗ | mesh->facets.delete_elements(to_delete); | |
| 552 | } | ||
| 553 | ✗ | } | |
| 554 | } | ||
| 555 | |||
| 556 | |||
| 557 | namespace GEO { | ||
| 558 | |||
| 559 | 4 | RVDCallback::RVDCallback() : | |
| 560 | 4 | seed_(index_t(-1)), | |
| 561 | 4 | simplex_(index_t(-1)), | |
| 562 | 4 | spinlocks_(nullptr) { | |
| 563 | 4 | } | |
| 564 | |||
| 565 | 8 | RVDCallback::~RVDCallback() { | |
| 566 | 8 | } | |
| 567 | |||
| 568 | ✗ | void RVDCallback::begin() { | |
| 569 | ✗ | seed_ = index_t(-1); | |
| 570 | ✗ | simplex_ = index_t(-1); | |
| 571 | ✗ | } | |
| 572 | |||
| 573 | ✗ | void RVDCallback::end() { | |
| 574 | ✗ | } | |
| 575 | |||
| 576 | /*********************************************************************/ | ||
| 577 | |||
| 578 | ✗ | RVDPolygonCallback::RVDPolygonCallback() { | |
| 579 | ✗ | } | |
| 580 | |||
| 581 | ✗ | RVDPolygonCallback::~RVDPolygonCallback() { | |
| 582 | ✗ | } | |
| 583 | |||
| 584 | ✗ | void RVDPolygonCallback::operator() ( | |
| 585 | index_t v, | ||
| 586 | index_t t, | ||
| 587 | const GEOGen::Polygon& C | ||
| 588 | ) const { | ||
| 589 | ✗ | const_cast<RVDPolygonCallback*>(this)->seed_ = v; | |
| 590 | ✗ | const_cast<RVDPolygonCallback*>(this)->simplex_ = t; | |
| 591 | ✗ | geo_argused(C); | |
| 592 | ✗ | } | |
| 593 | |||
| 594 | ✗ | void RVDPolygonCallback::begin() { | |
| 595 | ✗ | } | |
| 596 | |||
| 597 | ✗ | void RVDPolygonCallback::end() { | |
| 598 | ✗ | } | |
| 599 | |||
| 600 | /*********************************************************************/ | ||
| 601 | |||
| 602 | 4 | RVDPolyhedronCallback::RVDPolyhedronCallback() : | |
| 603 | 4 | facet_seed_(index_t(-1)), | |
| 604 | 4 | facet_tet_(index_t(-1)), | |
| 605 | 4 | last_seed_(index_t(-1)), | |
| 606 | 4 | simplify_internal_tet_facets_(false), | |
| 607 | 4 | simplify_voronoi_facets_(false), | |
| 608 | 4 | simplify_boundary_facets_(false), | |
| 609 | 4 | simplify_boundary_facets_angle_threshold_(0.0), | |
| 610 | 4 | tessellate_non_convex_facets_(false), | |
| 611 | 4 | use_mesh_(false), | |
| 612 | 4 | facet_is_skipped_(false), | |
| 613 |
4/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
|
4 | vertex_map_(nullptr) |
| 614 | { | ||
| 615 | 4 | } | |
| 616 | |||
| 617 | 8 | RVDPolyhedronCallback::~RVDPolyhedronCallback() { | |
| 618 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
|
8 | if(mesh_vertex_sym_.is_bound()) { |
| 619 | ✗ | mesh_vertex_sym_.unbind(); | |
| 620 | } | ||
| 621 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
|
8 | if(mesh_facet_seed_.is_bound()) { |
| 622 | ✗ | mesh_facet_seed_.unbind(); | |
| 623 | } | ||
| 624 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
|
8 | if(mesh_facet_tet_.is_bound()) { |
| 625 | ✗ | mesh_facet_tet_.unbind(); | |
| 626 | } | ||
| 627 | 8 | } | |
| 628 | |||
| 629 | ✗ | void RVDPolyhedronCallback::set_use_mesh(bool x) { | |
| 630 | ✗ | use_mesh_ = x; | |
| 631 | ✗ | if(!mesh_vertex_sym_.is_bound()) { | |
| 632 | ✗ | mesh_vertex_sym_.bind(mesh_.vertices.attributes(), "sym"); | |
| 633 | } | ||
| 634 | ✗ | if(!mesh_facet_seed_.is_bound()) { | |
| 635 | ✗ | mesh_facet_seed_.bind(mesh_.facets.attributes(),"seed"); | |
| 636 | } | ||
| 637 | ✗ | if(!mesh_facet_tet_.is_bound()) { | |
| 638 | ✗ | mesh_facet_tet_.bind(mesh_.facets.attributes(),"tet"); | |
| 639 | } | ||
| 640 | ✗ | } | |
| 641 | |||
| 642 | /********************************************************************/ | ||
| 643 | |||
| 644 | ✗ | void RVDPolyhedronCallback::begin_polyhedron( | |
| 645 | index_t seed, index_t tetrahedron | ||
| 646 | ) { | ||
| 647 | ✗ | geo_argused(seed); | |
| 648 | ✗ | geo_argused(tetrahedron); | |
| 649 | ✗ | } | |
| 650 | |||
| 651 | ✗ | void RVDPolyhedronCallback::begin_facet( | |
| 652 | index_t facet_seed, index_t facet_tet | ||
| 653 | ) { | ||
| 654 | ✗ | geo_argused(facet_seed); | |
| 655 | ✗ | geo_argused(facet_tet); | |
| 656 | ✗ | } | |
| 657 | |||
| 658 | ✗ | void RVDPolyhedronCallback::vertex( | |
| 659 | const double* geometry, const GEOGen::SymbolicVertex& symb | ||
| 660 | ) { | ||
| 661 | ✗ | geo_argused(geometry); | |
| 662 | ✗ | geo_argused(symb); | |
| 663 | ✗ | } | |
| 664 | |||
| 665 | ✗ | void RVDPolyhedronCallback::end_facet() { | |
| 666 | ✗ | } | |
| 667 | |||
| 668 | ✗ | void RVDPolyhedronCallback::end_polyhedron() { | |
| 669 | ✗ | } | |
| 670 | |||
| 671 | /********************************************************************/ | ||
| 672 | |||
| 673 | 1253 | void RVDPolyhedronCallback::begin_polyhedron_internal( | |
| 674 | index_t seed, index_t tetrahedron | ||
| 675 | ) { | ||
| 676 | 1253 | last_seed_ = seed; | |
| 677 | 1253 | seed_ = seed; | |
| 678 | 1253 | simplex_ = tetrahedron; | |
| 679 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1253 times.
|
1253 | if(use_mesh_) { |
| 680 | ✗ | vertex_map_ = new RVDVertexMap; | |
| 681 | } else { | ||
| 682 | 1253 | begin_polyhedron(seed, tetrahedron); | |
| 683 | } | ||
| 684 | 1253 | } | |
| 685 | |||
| 686 | 126029 | void RVDPolyhedronCallback::begin_facet_internal( | |
| 687 | index_t facet_seed, index_t facet_tet | ||
| 688 | ) { | ||
| 689 | 126029 | facet_seed_ = facet_seed; | |
| 690 | 126029 | facet_tet_ = facet_tet; | |
| 691 | 126029 | facet_is_skipped_ = ( | |
| 692 |
3/4✓ Branch 0 taken 126029 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 50032 times.
✓ Branch 3 taken 75997 times.
|
126029 | simplify_internal_tet_facets_ && (facet_tet != index_t(-1)) |
| 693 | ); | ||
| 694 |
2/2✓ Branch 0 taken 75997 times.
✓ Branch 1 taken 50032 times.
|
126029 | if(!facet_is_skipped_) { |
| 695 |
1/2✓ Branch 0 taken 75997 times.
✗ Branch 1 not taken.
|
75997 | if(use_mesh_) { |
| 696 | } else { | ||
| 697 | 75997 | begin_facet(facet_seed, facet_tet); | |
| 698 | } | ||
| 699 | } | ||
| 700 | 126029 | } | |
| 701 | |||
| 702 | 521718 | void RVDPolyhedronCallback::vertex_internal( | |
| 703 | const double* geometry, const GEOGen::SymbolicVertex& symb | ||
| 704 | ) { | ||
| 705 |
2/2✓ Branch 0 taken 302488 times.
✓ Branch 1 taken 219230 times.
|
521718 | if(!facet_is_skipped_) { |
| 706 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 302488 times.
|
302488 | if(use_mesh_) { |
| 707 | ✗ | index_t v = vertex_map_->find_or_create_vertex(seed(),symb); | |
| 708 | ✗ | if(v >= mesh_.vertices.nb()) { | |
| 709 | ✗ | mesh_.vertices.create_vertex(geometry); | |
| 710 | ✗ | mesh_vertex_sym_[v] = symb; | |
| 711 | } | ||
| 712 | ✗ | base_current_facet_.push_back(v); | |
| 713 | } else { | ||
| 714 | 302488 | vertex(geometry, symb); | |
| 715 | } | ||
| 716 | } | ||
| 717 | 521718 | } | |
| 718 | |||
| 719 | 126029 | void RVDPolyhedronCallback::end_facet_internal() { | |
| 720 |
2/2✓ Branch 0 taken 75997 times.
✓ Branch 1 taken 50032 times.
|
126029 | if(!facet_is_skipped_) { |
| 721 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 75997 times.
|
75997 | if(use_mesh_) { |
| 722 | ✗ | index_t f = mesh_.facets.nb(); | |
| 723 | ✗ | mesh_.facets.create_polygon(base_current_facet_.size()); | |
| 724 | ✗ | for(index_t i=0; i<base_current_facet_.size(); ++i) { | |
| 725 | ✗ | mesh_.facets.set_vertex(f, i, base_current_facet_[i]); | |
| 726 | } | ||
| 727 | ✗ | mesh_facet_seed_[f] = facet_seed(); | |
| 728 | ✗ | mesh_facet_tet_[f] = facet_tet(); | |
| 729 | ✗ | base_current_facet_.resize(0); | |
| 730 | } else { | ||
| 731 | 75997 | end_facet(); | |
| 732 | } | ||
| 733 | } | ||
| 734 | 126029 | facet_seed_ = index_t(-1); | |
| 735 | 126029 | facet_tet_ = index_t(-1); | |
| 736 | 126029 | } | |
| 737 | |||
| 738 | 1253 | void RVDPolyhedronCallback::end_polyhedron_internal() { | |
| 739 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1253 times.
|
1253 | if(use_mesh_) { |
| 740 | ✗ | mesh_.facets.connect(); | |
| 741 | ✗ | process_polyhedron_mesh(); | |
| 742 | ✗ | mesh_.clear(true,true); | |
| 743 | ✗ | delete vertex_map_; | |
| 744 | ✗ | vertex_map_ = nullptr; | |
| 745 | } else { | ||
| 746 | 1253 | end_polyhedron(); | |
| 747 | } | ||
| 748 | 1253 | seed_ = index_t(-1); | |
| 749 | 1253 | simplex_ = index_t(-1); | |
| 750 | 1253 | } | |
| 751 | |||
| 752 | /********************************************************************/ | ||
| 753 | |||
| 754 | ✗ | void RVDPolyhedronCallback::process_polyhedron_mesh() { | |
| 755 | ✗ | if(simplify_voronoi_facets_) { | |
| 756 | ✗ | simplify( | |
| 757 | ✗ | mesh_, | |
| 758 | ✗ | mesh_facet_seed_, | |
| 759 | simplify_boundary_facets_angle_threshold_ | ||
| 760 | ); | ||
| 761 | } | ||
| 762 | ✗ | if(tessellate_non_convex_facets_) { | |
| 763 | ✗ | tessellate_non_convex_facets(&mesh_); | |
| 764 | } | ||
| 765 | ✗ | begin_polyhedron(seed(), tet()); | |
| 766 | ✗ | for(index_t f=0; f<mesh_.facets.nb(); ++f) { | |
| 767 | ✗ | facet_seed_ = mesh_facet_seed_[f]; | |
| 768 | ✗ | facet_tet_ = mesh_facet_tet_[f]; | |
| 769 | ✗ | begin_facet(facet_seed_, facet_tet_); | |
| 770 | ✗ | for(index_t lv=0; lv<mesh_.facets.nb_vertices(f); ++lv) { | |
| 771 | ✗ | index_t v = mesh_.facets.vertex(f,lv); | |
| 772 | ✗ | vertex(mesh_.vertices.point_ptr(v), mesh_vertex_sym_[v]); | |
| 773 | } | ||
| 774 | ✗ | end_facet(); | |
| 775 | } | ||
| 776 | ✗ | end_polyhedron(); | |
| 777 | ✗ | } | |
| 778 | |||
| 779 | /********************************************************************/ | ||
| 780 | |||
| 781 | 4 | void RVDPolyhedronCallback::begin() { | |
| 782 | 4 | } | |
| 783 | |||
| 784 | 4 | void RVDPolyhedronCallback::end() { | |
| 785 | |||
| 786 | 4 | GEO::RVDPolyhedronCallback& callbacks = | |
| 787 | const_cast<GEO::RVDPolyhedronCallback&>(*this); | ||
| 788 | |||
| 789 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | if(simplify_internal_tet_facets_ && seed_ != index_t(-1)) { |
| 790 | 4 | callbacks.end_polyhedron_internal(); | |
| 791 | } | ||
| 792 | 4 | } | |
| 793 | |||
| 794 | 19547 | void RVDPolyhedronCallback::operator() ( | |
| 795 | index_t v, | ||
| 796 | index_t t, | ||
| 797 | const GEOGen::ConvexCell& C | ||
| 798 | ) const { | ||
| 799 | |||
| 800 | 19547 | GEO::RVDPolyhedronCallback& callbacks = | |
| 801 | const_cast<GEO::RVDPolyhedronCallback&>(*this); | ||
| 802 | |||
| 803 |
1/2✓ Branch 0 taken 19547 times.
✗ Branch 1 not taken.
|
19547 | if(simplify_internal_tet_facets_) { |
| 804 |
2/2✓ Branch 0 taken 1253 times.
✓ Branch 1 taken 18294 times.
|
19547 | if(v != last_seed_) { |
| 805 |
2/2✓ Branch 0 taken 1249 times.
✓ Branch 1 taken 4 times.
|
1253 | if(last_seed_ != index_t(-1)) { |
| 806 | 1249 | callbacks.end_polyhedron_internal(); | |
| 807 | } | ||
| 808 | 1253 | callbacks.begin_polyhedron_internal(v,t); | |
| 809 | } | ||
| 810 | } else { | ||
| 811 | ✗ | callbacks.begin_polyhedron_internal(v,t); | |
| 812 | } | ||
| 813 | |||
| 814 | // Remember that ConvexCell is represented in dual form ! | ||
| 815 | // - ConvexCell's vertices are facets | ||
| 816 | // - ConvexCell's triangles are vertices | ||
| 817 | |||
| 818 |
2/2✓ Branch 1 taken 309082 times.
✓ Branch 2 taken 19547 times.
|
328629 | for(index_t cv = 0; cv < C.max_v(); ++cv) { |
| 819 |
1/2✓ Branch 1 taken 309082 times.
✗ Branch 2 not taken.
|
309082 | signed_index_t ct = C.vertex_triangle(cv); |
| 820 |
2/2✓ Branch 0 taken 183053 times.
✓ Branch 1 taken 126029 times.
|
309082 | if(ct == -1) { |
| 821 | 183053 | continue; | |
| 822 | } | ||
| 823 |
2/8✓ Branch 1 taken 126029 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 126029 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
126029 | geo_debug_assert(C.triangle_is_used(index_t(ct))); |
| 824 | |||
| 825 |
1/2✓ Branch 1 taken 126029 times.
✗ Branch 2 not taken.
|
126029 | signed_index_t adjacent = C.vertex_id(cv); |
| 826 | 126029 | signed_index_t v_adj = -1; | |
| 827 | 126029 | signed_index_t t_adj = -1; | |
| 828 | |||
| 829 |
2/2✓ Branch 0 taken 50032 times.
✓ Branch 1 taken 75997 times.
|
126029 | if(adjacent < 0) { |
| 830 | // Negative adjacent indices correspond to | ||
| 831 | // tet-tet links | ||
| 832 | 50032 | t_adj = -adjacent - 1; | |
| 833 |
2/2✓ Branch 0 taken 71101 times.
✓ Branch 1 taken 4896 times.
|
75997 | } else if(adjacent > 0) { |
| 834 | // Positive adjacent indices correspond to | ||
| 835 | // Voronoi seed - Voronoi seed link | ||
| 836 | 71101 | v_adj = adjacent - 1; | |
| 837 | } // Zero adjacent indices corresponds to | ||
| 838 | // tet facet on border. | ||
| 839 | |||
| 840 |
1/2✓ Branch 1 taken 126029 times.
✗ Branch 2 not taken.
|
126029 | callbacks.begin_facet_internal(index_t(v_adj), index_t(t_adj)); |
| 841 | |||
| 842 | GEOGen::ConvexCell::Corner first( | ||
| 843 | index_t(ct), C.find_triangle_vertex(index_t(ct), cv) | ||
| 844 |
1/2✓ Branch 1 taken 126029 times.
✗ Branch 2 not taken.
|
126029 | ); |
| 845 | |||
| 846 | 126029 | GEOGen::ConvexCell::Corner c = first; | |
| 847 | do { | ||
| 848 |
1/2✓ Branch 1 taken 521718 times.
✗ Branch 2 not taken.
|
521718 | const GEOGen::Vertex& vx = C.triangle_dual(c.t); |
| 849 |
1/2✓ Branch 3 taken 521718 times.
✗ Branch 4 not taken.
|
521718 | callbacks.vertex_internal(vx.point(), vx.sym()); |
| 850 |
1/2✓ Branch 1 taken 521718 times.
✗ Branch 2 not taken.
|
521718 | C.move_to_next_around_vertex(c); |
| 851 |
2/2✓ Branch 1 taken 395689 times.
✓ Branch 2 taken 126029 times.
|
521718 | } while(c != first); |
| 852 |
1/2✓ Branch 1 taken 126029 times.
✗ Branch 2 not taken.
|
126029 | callbacks.end_facet_internal(); |
| 853 | } | ||
| 854 | |||
| 855 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 19547 times.
|
19547 | if(!simplify_internal_tet_facets_) { |
| 856 | ✗ | callbacks.end_polyhedron_internal(); | |
| 857 | } | ||
| 858 | 19547 | } | |
| 859 | |||
| 860 | /********************************************************************/ | ||
| 861 | |||
| 862 | 4 | BuildRVDMesh::BuildRVDMesh(Mesh& output_mesh) : | |
| 863 |
4/8✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
|
4 | output_mesh_(output_mesh), shrink_(0.0) { |
| 864 | 4 | cell_vertex_map_ = nullptr; | |
| 865 | 4 | global_vertex_map_ = nullptr; | |
| 866 | 4 | current_cell_id_ = 0; | |
| 867 | 4 | generate_ids_ = false; | |
| 868 | 4 | } | |
| 869 | |||
| 870 | 8 | BuildRVDMesh::~BuildRVDMesh() { | |
| 871 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
8 | if(generate_ids_) { |
| 872 | ✗ | cell_id_.unbind(); | |
| 873 | ✗ | seed_id_.unbind(); | |
| 874 | ✗ | vertex_id_.unbind(); | |
| 875 | ✗ | facet_seed_id_.unbind(); | |
| 876 | ✗ | delete global_vertex_map_; | |
| 877 | ✗ | global_vertex_map_ = nullptr; | |
| 878 | } | ||
| 879 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
8 | delete cell_vertex_map_; |
| 880 | 8 | cell_vertex_map_ = nullptr; | |
| 881 | 8 | } | |
| 882 | |||
| 883 | 4 | void BuildRVDMesh::set_generate_ids(bool x) { | |
| 884 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if(x == generate_ids_) { |
| 885 | 4 | return; | |
| 886 | } | ||
| 887 | ✗ | generate_ids_ = x; | |
| 888 | ✗ | if(generate_ids_) { | |
| 889 | ✗ | cell_id_.bind( | |
| 890 | ✗ | output_mesh_.facets.attributes(), "cell_id" | |
| 891 | ); | ||
| 892 | ✗ | seed_id_.bind( | |
| 893 | ✗ | output_mesh_.facets.attributes(), "seed_id" | |
| 894 | ); | ||
| 895 | ✗ | vertex_id_.bind( | |
| 896 | ✗ | output_mesh_.vertices.attributes(), "vertex_id" | |
| 897 | ); | ||
| 898 | ✗ | facet_seed_id_.bind( | |
| 899 | ✗ | output_mesh_.facets.attributes(), "facet_seed_id" | |
| 900 | ); | ||
| 901 | ✗ | global_vertex_map_ = new RVDVertexMap; | |
| 902 | } else { | ||
| 903 | ✗ | cell_id_.unbind(); | |
| 904 | ✗ | seed_id_.unbind(); | |
| 905 | ✗ | vertex_id_.unbind(); | |
| 906 | ✗ | facet_seed_id_.unbind(); | |
| 907 | ✗ | delete global_vertex_map_; | |
| 908 | ✗ | global_vertex_map_ = nullptr; | |
| 909 | } | ||
| 910 | } | ||
| 911 | |||
| 912 | 4 | void BuildRVDMesh::set_shrink(double x) { | |
| 913 | 4 | shrink_ = x; | |
| 914 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if(shrink_ != 0.0) { |
| 915 | ✗ | set_use_mesh(true); | |
| 916 | } | ||
| 917 | 4 | } | |
| 918 | |||
| 919 | 4 | void BuildRVDMesh::begin() { | |
| 920 | 4 | RVDPolyhedronCallback::begin(); | |
| 921 | 4 | output_mesh_.clear(); | |
| 922 | 4 | output_mesh_.vertices.set_dimension(3); | |
| 923 | 4 | } | |
| 924 | |||
| 925 | 4 | void BuildRVDMesh::end() { | |
| 926 | 4 | RVDPolyhedronCallback::end(); | |
| 927 | 4 | output_mesh_.facets.connect(); | |
| 928 | 4 | } | |
| 929 | |||
| 930 | |||
| 931 | 1253 | void BuildRVDMesh::begin_polyhedron(index_t seed, index_t tetrahedron) { | |
| 932 | 1253 | geo_argused(tetrahedron); | |
| 933 | 1253 | geo_argused(seed); | |
| 934 |
2/2✓ Branch 0 taken 1249 times.
✓ Branch 1 taken 4 times.
|
1253 | delete cell_vertex_map_; |
| 935 |
2/6✓ Branch 2 taken 1253 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1253 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1253 | cell_vertex_map_ = new RVDVertexMap; |
| 936 | 1253 | cell_vertex_map_->set_first_vertex_index( | |
| 937 | 1253 | output_mesh_.vertices.nb() | |
| 938 | ); | ||
| 939 | 1253 | } | |
| 940 | |||
| 941 | 75997 | void BuildRVDMesh::begin_facet(index_t facet_seed, index_t facet_tet_facet) { | |
| 942 | 75997 | geo_argused(facet_seed); | |
| 943 | 75997 | geo_argused(facet_tet_facet); | |
| 944 | 75997 | current_facet_.resize(0); | |
| 945 | 75997 | } | |
| 946 | |||
| 947 | 302488 | void BuildRVDMesh::vertex( | |
| 948 | const double* geometry, const GEOGen::SymbolicVertex& symb | ||
| 949 | ) { | ||
| 950 |
1/2✓ Branch 2 taken 302488 times.
✗ Branch 3 not taken.
|
302488 | index_t v = cell_vertex_map_->find_or_create_vertex(seed(), symb); |
| 951 |
2/2✓ Branch 1 taken 77410 times.
✓ Branch 2 taken 225078 times.
|
302488 | if(v >= output_mesh_.vertices.nb()) { |
| 952 |
1/2✓ Branch 1 taken 77410 times.
✗ Branch 2 not taken.
|
77410 | output_mesh_.vertices.create_vertex(geometry); |
| 953 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 77410 times.
|
77410 | if(generate_ids_) { |
| 954 | ✗ | vertex_id_[v] = int( | |
| 955 | ✗ | global_vertex_map_->find_or_create_vertex(seed(), symb) | |
| 956 | ); | ||
| 957 | } | ||
| 958 | } | ||
| 959 |
1/2✓ Branch 1 taken 302488 times.
✗ Branch 2 not taken.
|
302488 | current_facet_.push_back(v); |
| 960 | 302488 | } | |
| 961 | |||
| 962 | 75997 | void BuildRVDMesh::end_facet() { | |
| 963 | 75997 | index_t f = output_mesh_.facets.nb(); | |
| 964 | 75997 | output_mesh_.facets.create_polygon(current_facet_.size()); | |
| 965 |
2/2✓ Branch 1 taken 302488 times.
✓ Branch 2 taken 75997 times.
|
378485 | for(index_t i=0; i<current_facet_.size(); ++i) { |
| 966 | 302488 | output_mesh_.facets.set_vertex(f,i,current_facet_[i]); | |
| 967 | } | ||
| 968 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 75997 times.
|
75997 | if(generate_ids_) { |
| 969 | ✗ | seed_id_[f] = int(seed()); | |
| 970 | ✗ | cell_id_[f] = int(current_cell_id_); | |
| 971 | ✗ | facet_seed_id_[f] = int(facet_seed()); | |
| 972 | } | ||
| 973 | 75997 | } | |
| 974 | |||
| 975 | 1253 | void BuildRVDMesh::end_polyhedron() { | |
| 976 | 1253 | ++current_cell_id_; | |
| 977 | 1253 | } | |
| 978 | |||
| 979 | ✗ | void BuildRVDMesh::process_polyhedron_mesh() { | |
| 980 | ✗ | if(shrink_ != 0.0 && mesh_.vertices.nb() != 0) { | |
| 981 | ✗ | vec3 center(0.0, 0.0, 0.0); | |
| 982 | ✗ | for(index_t v=0; v<mesh_.vertices.nb(); ++v) { | |
| 983 | ✗ | center += vec3(mesh_.vertices.point_ptr(v)); | |
| 984 | } | ||
| 985 | ✗ | center = (1.0 / double(mesh_.vertices.nb())) * center; | |
| 986 | ✗ | for(index_t v=0; v<mesh_.vertices.nb(); ++v) { | |
| 987 | ✗ | vec3 p(mesh_.vertices.point_ptr(v)); | |
| 988 | ✗ | p = shrink_ * center + (1.0 - shrink_) * p; | |
| 989 | ✗ | mesh_.vertices.point_ptr(v)[0] = p.x; | |
| 990 | ✗ | mesh_.vertices.point_ptr(v)[1] = p.y; | |
| 991 | ✗ | mesh_.vertices.point_ptr(v)[2] = p.z; | |
| 992 | } | ||
| 993 | } | ||
| 994 | ✗ | RVDPolyhedronCallback::process_polyhedron_mesh(); | |
| 995 | ✗ | } | |
| 996 | |||
| 997 | } | ||
| 998 |