1 2 #include <exploragram/hexdom/quad_dominant.h> 3 #include <exploragram/hexdom/meshcomesh.h> 4 #include <exploragram/hexdom/PGP.h> 5 #include <exploragram/hexdom/basic.h> 6 #include <exploragram/hexdom/extra_connectivity.h> 7 #include <exploragram/hexdom/mesh_inspector.h> 8 #include <exploragram/hexdom/intersect_tools.h> 9 #include <exploragram/hexdom/polygon.h> 10 11 #include <geogram/numerics/matrix_util.h> 12 #include <geogram/basic/permutation.h> 13 #include <geogram/mesh/triangle_intersection.h> 14 #include <geogram/mesh/mesh_tetrahedralize.h> 15 #include <geogram/delaunay/delaunay.h> 16 #include <geogram/points/nn_search.h> 17 #include <geogram/points/colocate.h> 18 #include <geogram/mesh/mesh_repair.h> 19 #include <geogram/mesh/mesh_fill_holes.h> 20 #include <geogram/mesh/mesh_geometry.h> 21 22 23 #include <algorithm> 24 #include <queue> 25 #include <stack> 26 #include <map> 27 28 namespace GEO { 29 export_boundary_with_uv(Mesh * m,Mesh * hex,const char * uv_name,const char * singular_name)30 void export_boundary_with_uv(Mesh* m, Mesh* hex, const char* uv_name, const char* singular_name) { 31 Attribute<GEO::vec2> uv(hex->facet_corners.attributes(), uv_name); 32 Attribute<index_t> singtri(hex->facets.attributes(), singular_name); 33 34 Attribute<bool> has_param(m->cell_facets.attributes(), "has_param"); 35 Attribute<vec3> UC(m->cell_corners.attributes(), "U"); 36 37 // STEP 1: compute number of boundary vertices AND create a lookup table tet mesh vertex -> boundary surface vertex 38 vector<index_t> tetV_to_facetV(m->vertices.nb(), NOT_AN_ID); 39 index_t nb_boundary_V = 0; 40 FOR(c, m->cells.nb()) FOR(cf, m->cells.nb_facets(c)) { 41 if (NO_CELL != m->cells.adjacent(c, cf)) continue; 42 FOR(cfv, m->cells.facet_nb_vertices(c, cf)) { 43 index_t v = m->cells.facet_vertex(c, cf, cfv); 44 if (NOT_AN_ID == tetV_to_facetV[v]) { 45 tetV_to_facetV[v] = nb_boundary_V++; 46 } 47 } 48 } 49 50 // STEP 2: copy boundary vertices to the new mesh 51 hex->vertices.create_vertices(nb_boundary_V); 52 FOR(v, m->vertices.nb()) 53 if (NOT_AN_ID != tetV_to_facetV[v]) 54 hex->vertices.point(tetV_to_facetV[v]) = m->vertices.point(v); 55 56 // STEP 3: extract triangles and their parameterization (when possible) 57 FOR(c, m->cells.nb()) FOR(cf, m->cells.nb_facets(c)) { 58 if (m->cells.adjacent(c, cf) != NO_CELL) continue; 59 index_t tet_verts[3]; 60 FOR(cfv, 3) { 61 tet_verts[cfv] = m->cells.facet_vertex(c, cf, cfv); 62 } 63 64 index_t tet_corners[3]; 65 FOR(cfv, 3) FOR(test, 4) if (m->cell_corners.vertex(m->cells.corner(c, test)) == tet_verts[cfv]) 66 tet_corners[cfv] = m->cells.corner(c, test); 67 68 69 vec3 lX[3], lU[3]; 70 FOR(cfv, 3) { 71 lX[cfv] = m->vertices.point(tet_verts[cfv]); 72 lU[cfv] = UC[tet_corners[cfv]]; 73 } 74 75 // STEP 2.2: non singular boundary triangles are easy to extract 76 index_t fid = hex->facets.create_triangle( 77 tetV_to_facetV[tet_verts[0]], 78 tetV_to_facetV[tet_verts[1]], 79 tetV_to_facetV[tet_verts[2]] 80 ); 81 82 bool has_valid_2d_param = false; 83 if (has_param[m->cells.facet(c, cf)]) { 84 FOR(dim, 3) { // we are looking for the param dimension that goes inside the volume 85 if (lU[0][dim] != lU[1][dim] || lU[0][dim] != lU[2][dim]) continue; 86 has_valid_2d_param = true; 87 FOR(lv, 3) { 88 uv[hex->facets.corner(fid, lv)] = vec2(lU[lv][(dim + 1) % 3], lU[lv][(dim + 2) % 3]); 89 } 90 } 91 #if 0 92 // mirror the copy when necessary 93 if (det(uv[hex->facets.corner(fid, 1)] - uv[hex->facets.corner(fid, 0)], 94 uv[hex->facets.corner(fid, 2)] - uv[hex->facets.corner(fid, 0)]) < 0) 95 FOR(lv, 3) uv[hex->facets.corner(fid, lv)][0] *= -1.; 96 #endif 97 } 98 99 if (has_valid_2d_param) { // check param quality and invalidate it, if necessary 100 TrglGradient grd(lX[0], lX[1], lX[2]); 101 vec3 grduv[2]; 102 FOR(d, 2) grduv[d] = grd.gradient_3d(uv[hex->facets.corner(fid, 0)][d], uv[hex->facets.corner(fid, 1)][d], uv[hex->facets.corner(fid, 2)][d]); 103 FOR(d, 2) FOR(dd, 3) if (Numeric::is_nan(grduv[d][dd])) has_valid_2d_param = false; 104 #if 0 105 if (grduv[0].length() > 10. * grduv[1].length()) has_valid_2d_param = false; 106 if (grduv[1].length() > 10. * grduv[0].length()) has_valid_2d_param = false; 107 if (std::abs(dot(normalize(grduv[0]), normalize(grduv[1]))) > cos(M_PI / 4.)) has_valid_2d_param = false; 108 #endif 109 } 110 111 singtri[fid] = !has_valid_2d_param; 112 } 113 } 114 triangulate_surface_preserve_attributes(Mesh * m)115 static bool triangulate_surface_preserve_attributes(Mesh* m) { 116 bool result = true; 117 plop("triangulating surface after embedding isoUV"); 118 vector<index_t> to_kill(m->facets.nb(), false); 119 FOR(f, m->facets.nb()) { 120 geo_assert(m->facets.nb_corners(f) >= 3); 121 if (m->facets.nb_corners(f) == 3) continue; 122 123 vector<vec3> pts; 124 FOR(lc, m->facets.nb_corners(f)) { 125 pts.push_back(X(m)[m->facets.vertex(f, lc)]); 126 } 127 128 vector<index_t> triangles; 129 bool success = Poly3d(pts).try_triangulate_minweight(triangles); 130 result = result && success; 131 132 geo_assert(0 == triangles.size() % 3); 133 134 FOR(new_face, triangles.size() / 3) { 135 index_t new_f = m->facets.create_triangle( 136 m->facets.vertex(f, triangles[3 * new_face + 0]), 137 m->facets.vertex(f, triangles[3 * new_face + 1]), 138 m->facets.vertex(f, triangles[3 * new_face + 2])); 139 140 to_kill.push_back(false); 141 m->facets.attributes().copy_item(new_f, f); 142 FOR(nfc,3) { 143 m->facet_corners.attributes().copy_item(m->facets.corner(new_f, nfc), m->facets.corner(f, triangles[3 * new_face + nfc])); 144 } 145 } 146 to_kill[f] = true; 147 } 148 m->facets.delete_elements(to_kill); 149 return result; 150 } 151 152 153 /** 154 * INPUT: facets with xyz geometry and uv coordinates s.t. no edge is of size 0 155 singular bool per triangle 156 * integer values of uv --- interpolated along edge 'e' --- matches on both side of 'e' 157 on presuppose que les U ont �t� pre-snapp�s sur la grille entiere 158 * OUTPUT: facets with uv coordinates s.t. 159 * a new vertex is inserted (with interpolated xyz and uv) at every integer values of uv along edges 160 * nothing prevents some vertices around a facet to share the same uv's 161 ATTENTION: it modifies the parameterization when snaps grid corners on edges! 162 */ 163 split_edges_by_iso_uvs(Mesh * m,const char * uv_name,const char * singular_name)164 void split_edges_by_iso_uvs(Mesh* m, const char *uv_name, const char *singular_name) { 165 Attribute<GEO::vec2> uv(m->facet_corners.attributes(), uv_name); 166 Attribute<index_t> singular(m->facets.attributes(), singular_name); 167 168 index_t nb_init_facets = m->facets.nb(); 169 Attribute<bool> added_vertices(m->vertices.attributes(), "added_vertices"); 170 Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid"); 171 Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet"); 172 173 typedef std::pair<index_t /*v_id*/, vec2 /*uv*/> NewCorner; 174 vector<vector<NewCorner> > new_corners(m->facet_corners.nb()); 175 176 reach("create vertices and init arrays of NewCorner to insert on edges"); 177 { 178 FacetsExtraConnectivity fec(m); 179 FOR(h, m->facet_corners.nb()) { 180 index_t opp = fec.opposite(h); 181 182 if (singular[fec.facet(h)]) continue; // the face is responsible for the edge if opp<h or if its opposite is singular 183 if (NOT_AN_ID != opp && !singular[fec.facet(opp)] && opp > h) continue; 184 185 vec2 lU[2] = { uv[h], uv[fec.next(h)] }; 186 187 #if 0 188 bool trans_func_found = false; 189 index_t R; 190 vec2 T; 191 if (NOT_AN_ID != opp && !singular[fec.facet(opp)]) { 192 vec2 A = lU[1]-lU[0], B = uv[opp]-uv[fec.next(opp)], C = lU[0]; 193 if (std::abs(A.length()-B.length())<1e-3) { 194 double mindist = std::numeric_limits<double>::max(); 195 FOR(r,4) { 196 if ((A-B).length()<mindist) { 197 R = r; 198 mindist = (A-B).length(); 199 } 200 A = vec2(-A[1], A[0]); 201 } 202 FOR(r, R) { 203 A = vec2(-A[1], A[0]); 204 C = vec2(-C[1], C[0]); 205 } 206 plop(mindist); 207 vec2 Ttmp = uv[fec.next(opp)] - C; 208 T = vec2(round(Ttmp[0]), round(Ttmp[1])); 209 trans_func_found = (mindist<1e-3 && (T-Ttmp).length()<1e-3); 210 } else { 211 plop("GNA?!"); 212 } 213 if (!trans_func_found && NOT_AN_ID!=opp && !singular[fec.facet(opp)]) { 214 plop("ATTENTION, grids do not match"); 215 singular[fec.facet(h)] = true; 216 singular[fec.facet(opp)] = true; 217 continue; 218 } 219 } 220 #endif 221 vector<double> coeff; 222 FOR(coord, 2) { // find barycentric coordinates of both u and v integer isos 223 double v[2] = { lU[0][coord], lU[1][coord] }; 224 double from = floor(std::min(v[0], v[1])) + 1; 225 double to = std::max(v[0], v[1]); 226 if (to-from > 1000) continue; 227 for (double iso = from; iso < to; iso += 1.) { 228 double c = (iso - v[0]) / (v[1] - v[0]); // v[0] is far from v[1] (U was pre-snapped to integers with .05 tolerance) 229 if (!Numeric::is_nan(c) && c > 0 && c < 1) { 230 // vec2 u = lU[0] + c*(lU[1] - lU[0]); 231 // u[coord] = iso; 232 coeff.push_back(c); 233 } 234 } 235 } 236 std::sort(coeff.begin(), coeff.end(), std::less<double>()); // we need to sort in order to merge close values 237 238 vector<vec3> pts; 239 vector<vec2> lu; 240 vector<vec2> luopp; 241 242 FOR(i, coeff.size()) { 243 // a + c*(b-a) returns exactly a when a==b and no NaNs involved 244 // no need to worry about the cases when c==0 and c==1, because of the presnapping of the parameterization 245 vec3 pt = X(m)[fec.org(h)] + coeff[i] * (X(m)[fec.dest(h)] - X(m)[fec.org(h)]); 246 vec2 u = lU[0] + coeff[i] * (lU[1] - lU[0]); 247 vec2 uopp; 248 if (NOT_AN_ID != opp) uopp = uv[fec.next(opp)] + coeff[i] * (uv[opp] - uv[fec.next(opp)]); 249 250 FOR(d, 2) { // we must guarantee that new vertices have (at least) one integer component 251 if (std::abs(u[d] - round(u[d])) < 1e-10) { 252 u[d] = round(u[d]); 253 } 254 255 if (NOT_AN_ID != opp && std::abs(uopp[d] - round(uopp[d])) < 1e-10) { 256 uopp[d] = round(uopp[d]); 257 } 258 } 259 260 pts.push_back(pt); 261 lu.push_back(u); 262 if (NOT_AN_ID != opp) luopp.push_back(uopp); 263 } 264 265 // create vertices 266 index_t off = m->vertices.create_vertices(pts.size()); 267 FOR(i, pts.size()) { 268 added_vertices[i+off] = true; 269 resp_facet[i+off] = orig_tri_fid[fec.facet(h)]; 270 } 271 FOR(i, pts.size()) { 272 m->vertices.point(off + i) = pts[i]; 273 new_corners[h].push_back(NewCorner(off + i, lu[i])); 274 if (NOT_AN_ID != opp) 275 new_corners[opp].push_back(NewCorner(off + i, luopp[i])); 276 } 277 if (NOT_AN_ID != opp) 278 std::reverse(new_corners[opp].begin(), new_corners[opp].end()); 279 } 280 } 281 282 // now we have all the corners to insert, we create new facets for all the surface, old facets are to delete 283 284 reach("split edges"); 285 FOR(f, nb_init_facets) { 286 vector<index_t> polyV; 287 vector<vec2> poly_uv; 288 FOR(fc, m->facets.nb_corners(f)) { 289 polyV.push_back(m->facets.vertex(f, fc)); 290 index_t c = m->facets.corner(f, fc); 291 poly_uv.push_back(uv[c]); 292 FOR(i, new_corners[c].size()) { 293 polyV.push_back(new_corners[c][i].first); 294 poly_uv.push_back(new_corners[c][i].second); 295 } 296 } 297 index_t nf = m->facets.create_polygon(polyV); 298 m->facets.attributes().copy_item(nf ,f); 299 FOR(fc, m->facets.nb_corners(nf)) 300 uv[m->facets.corner(nf, fc)] = poly_uv[fc]; 301 } 302 303 reach("kill facets"); 304 vector<index_t> to_kill(nb_init_facets, true); // kill old (pre-split) facets 305 to_kill.resize(m->facets.nb(), false); 306 m->facets.delete_elements(to_kill); 307 } 308 309 310 /** 311 * INPUT: facets with uv coordinates, where many vertices have integer values in uv 312 * OUTPUT: facets with uv coordinates, inclusing edges that have one coordinate of uv that is constant and integer valued 313 * 314 * remark: some polylines of these new edges are likely to be pre-images of edges of the regular grid... 315 */ 316 317 find_degenerate_facets(Mesh * m,vector<index_t> & degenerate)318 void find_degenerate_facets(Mesh* m, vector<index_t> °enerate) { 319 degenerate.clear(); 320 FOR(f, m->facets.nb()) { 321 index_t nbv = m->facets.nb_corners(f); 322 // geo_assert(3 == m->facets.nb_corners(f)); 323 FOR (c1, nbv) { 324 FOR (c2, nbv) { 325 if (c1 == c2) continue; 326 index_t v1 = m->facets.vertex(f, c1); 327 index_t v2 = m->facets.vertex(f, c2); 328 329 geo_assert( v1!= v2 ); 330 if (X(m)[v1][0] == X(m)[v2][0] && X(m)[v1][1] == X(m)[v2][1] && X(m)[v1][2] == X(m)[v2][2]) { 331 degenerate.push_back(f); 332 } 333 } 334 } 335 } 336 } 337 imprint(Mesh * m,const char * uv_name,const char * singular_name)338 void imprint(Mesh* m, const char *uv_name, const char *singular_name) { 339 { 340 Attribute<index_t> singular(m->facets.attributes(), singular_name); 341 Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet"); 342 Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid"); 343 FOR(i, m->facets.nb() ) orig_tri_fid[i] = i; 344 FOR(v, m->vertices.nb()) resp_facet[v] = NOT_AN_ID; 345 346 check_no_intersecting_faces(m); 347 } 348 349 350 351 Mesh m_bak; 352 m_bak.copy(*m); 353 354 for(;;) { 355 split_edges_by_iso_uvs(m, uv_name, singular_name); 356 facets_split(m, uv_name, singular_name); 357 358 triangulate_surface_preserve_attributes(m); 359 360 vector<index_t> fails; 361 find_degenerate_facets(m, fails); 362 363 vector<index_t> intersections; 364 find_self_intersections(m, intersections); 365 fails.insert(fails.end(), intersections.begin(), intersections.end()); // degenerate triangles or intersecting ones, all the same, I do not want them 366 if (!fails.size()) break; 367 368 plop("imprint fail, need to undo"); 369 { 370 Attribute<index_t> singular(m_bak.facets.attributes(), singular_name); 371 Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet"); 372 373 FOR(i, fails.size()) { 374 FOR(j, 3) { 375 index_t f = resp_facet[m->facets.vertex(fails[i], j)]; 376 if (NOT_AN_ID!=f) singular[f] = true; 377 } 378 } 379 } 380 m->copy(m_bak); 381 } 382 } 383 facets_split(Mesh * m,const char * uv_name,const char * singular_name)384 void facets_split(Mesh* m, const char *uv_name, const char *singular_name) { 385 Attribute<bool> added_vertices(m->vertices.attributes(), "added_vertices"); 386 Attribute<GEO::vec2> uv(m->facet_corners.attributes(), uv_name); 387 Attribute<index_t> singular(m->facets.attributes(), singular_name); 388 389 Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid"); 390 Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet"); 391 392 vector<index_t> to_kill(m->facets.nb(), 0); 393 FOR(f, m->facets.nb()) { 394 if (singular[f]) continue; 395 index_t nbc = m->facets.nb_corners(f); 396 397 398 // STEP 1: find a couple (coordinate, iso value) to split the facet 399 index_t coord = index_t(-1); 400 double iso = 0; 401 FOR(test_coord, 2) { 402 double min_v = 1e20; 403 double max_v = -1e20; 404 FOR(fc, nbc) { 405 min_v = std::min(min_v, uv[m->facets.corner(f, fc)][test_coord]); 406 max_v = std::max(max_v, uv[m->facets.corner(f, fc)][test_coord]); 407 } 408 if (floor(min_v)+1 < max_v) { // floor(min_v)+1 is the first integer value strictly superior to min_v 409 coord = test_coord; 410 iso = floor(min_v)+1; 411 break; 412 } 413 } 414 if (coord == index_t(-1)) continue; 415 416 // STEP 2: if STEP 1 succedeed, compute the extremities of the new edge 417 index_t cut[2] = { NOT_AN_ID, NOT_AN_ID }; 418 419 FOR(fc, nbc) { 420 if (uv[m->facets.corner(f, fc)][coord] != iso) continue; 421 if (cut[0] == NOT_AN_ID) { 422 cut[0] = fc; 423 } else if (fc != next_mod(cut[0], nbc) && next_mod(fc, nbc) != cut[0]) { // no biangles 424 cut[1] = fc; 425 } 426 } 427 if (cut[1] == NOT_AN_ID) continue; 428 429 { // let us check that the cut separates the facet 430 bool inf1=true, sup1=true, inf2=true, sup2=true; 431 for (index_t fc = cut[0]+1; fc<cut[1]; fc++) { 432 double tex = uv[m->facets.corner(f, fc)][coord]; 433 inf1 = inf1 && (tex<iso); 434 sup1 = sup1 && (tex>iso); 435 } 436 for (index_t fc = cut[1]+1; fc<cut[0]+nbc; fc++) { 437 double tex = uv[m->facets.corner(f, fc%nbc)][coord]; 438 inf2 = inf2 && (tex<iso); 439 sup2 = sup2 && (tex>iso); 440 } 441 if (!( (inf1 && sup2) || (sup1 && inf2) )) { 442 plop("WARNING: facet cannot be properly cut by the iso"); 443 continue; 444 } 445 } 446 447 448 to_kill[f] = true; 449 450 // STEP 3: compute new vertices (iso-integer value of uv) on the new edge 451 vector<vec3> nv_pts; 452 vector<vec2> nv_uv; 453 { 454 vec3 lX[2]; 455 FOR(i, 2) lX[i] = m->vertices.point(m->facets.vertex(f, cut[i])); 456 vec2 lU[2]; 457 FOR(i, 2) lU[i] = uv[m->facets.corner(f, cut[i])]; 458 vector<double> coeff; 459 double v[2] = { lU[0][(coord+1)%2], lU[1][(coord+1)%2] }; // recall that coord is the cutting dimension 460 461 for (double cur_iso = ceil(std::min(v[0], v[1])); cur_iso < std::max(v[0], v[1]); cur_iso += 1.0) { 462 double c = (cur_iso - v[0]) / (v[1] - v[0]); 463 if (!Numeric::is_nan(c) && c > 0 && c < 1) 464 coeff.push_back(c); 465 } 466 467 std::sort(coeff.begin(), coeff.end(), std::less<double>()); 468 FOR(i, coeff.size()) { 469 vec3 x = lX[0] + coeff[i] * (lX[1] - lX[0]); // it guarantees x==lX[0] when lX[0]==lX[1] 470 vec2 u = lU[0] + coeff[i] * (lU[1] - lU[0]); // no need to worry about coeff[i]==0 and coeff[i]==1 because of the parameterization pre-snapping 471 nv_pts.push_back(x); 472 nv_uv.push_back(u); 473 } 474 // new vertices must have only integer values of uv --- remove possible numerical imprecision 475 FOR(i, nv_pts.size()) { 476 FOR(d, 2) { 477 nv_uv[i][d] = round(nv_uv[i][d]); 478 } 479 } 480 } 481 482 // STEP 4: create new vertices and new faces 483 index_t off = m->vertices.create_vertices(nv_pts.size()); 484 FOR(i, nv_pts.size()) { 485 resp_facet[off+i] = orig_tri_fid[f]; 486 added_vertices[off+i] = true; 487 X(m)[off + i] = nv_pts[i]; 488 } 489 FOR(half, 2) { 490 vector <index_t> lv; 491 vector <vec2> luv; 492 493 // add original vertices 494 index_t cir = cut[half]; 495 do { 496 lv.push_back(m->facets.vertex(f, cir)); 497 luv.push_back(uv[m->facets.corner(f, cir)]); 498 cir = next_mod(cir, nbc); 499 } while (cir != cut[(half + 1) % 2]); 500 lv.push_back(m->facets.vertex(f, cir)); 501 luv.push_back(uv[m->facets.corner(f, cir)]); 502 503 // add new vertices 504 FOR(i, nv_pts.size()) { 505 index_t ind = i; 506 if (half == 0) ind = nv_pts.size() - 1 - i; 507 lv.push_back(off + ind); 508 luv.push_back(nv_uv[ind]); 509 } 510 511 index_t fid = m->facets.create_polygon(lv); 512 m->facets.attributes().copy_item(fid, f); 513 FOR(fc, m->facets.nb_corners(fid)) { 514 uv[m->facets.corner(fid, fc)] = luv[fc]; 515 } 516 to_kill.push_back(false); 517 } 518 } 519 m->facets.delete_elements(to_kill); 520 } 521 522 /** 523 * INPUT: facets with uv coordinates 524 * OUTPUT: chart facet attribute s.t. the chart frontier is included in edges that are iso-integer value of 'uv' 525 */ 526 indir_root(index_t i,vector<index_t> & indir)527 inline index_t indir_root(index_t i, vector<index_t>& indir) { 528 while (i != indir[i]) i = indir[i]; 529 return i; 530 } 531 532 // merges charts for adjacent facets under two conditions: 533 // 1) both facets are not singular 534 // 2) the shared edge is not integer iso in both facets mark_charts(Mesh * m,const char * uv_name,const char * chart_name,const char * singular_name)535 void mark_charts(Mesh* m, const char *uv_name, const char *chart_name, const char *singular_name) { // 2-manifold surface is supposed 536 Attribute<bool> isovalue(m->facet_corners.attributes(), "isovalue"); 537 Attribute<bool> quadelement(m->facets.attributes(), "quadelement"); 538 Attribute<bool> quadcorners(m->vertices.attributes(), "quadcorners"); 539 Attribute<bool> added_vertices(m->vertices.attributes(), "added_vertices"); 540 541 Attribute<GEO::vec2> uv(m->facet_corners.attributes(), uv_name); 542 Attribute<index_t> singular(m->facets.attributes(), singular_name); 543 Attribute<index_t> chart(m->facets.attributes(), chart_name); 544 545 FacetsExtraConnectivity fec(m); 546 547 vector<index_t> indir(m->facets.nb()); 548 { // fill isovalue edge attribute and merge quad candidate charts 549 FOR(f, m->facets.nb()) { 550 indir[f] = f; 551 } 552 553 FOR(h, m->facet_corners.nb()) { 554 index_t hopp = fec.opposite(h); 555 index_t f = fec.facet(h); 556 index_t fopp = NOT_AN_ID==hopp ? NOT_AN_ID : fec.facet(hopp); 557 558 if (singular[f] || (NOT_AN_ID!=hopp && hopp>h)) continue; 559 560 bool cut = false; 561 index_t test[2] = { h, hopp }; 562 FOR(lh, 2) { 563 if (lh && (NOT_AN_ID==hopp || singular[fopp])) break; 564 FOR(coord, 2) { 565 cut = cut || (uv[test[lh]][coord] == uv[fec.next(test[lh])][coord] && uv[test[lh]][coord] == round(uv[test[lh]][coord])); 566 } 567 } 568 569 if (cut) { 570 isovalue[h] = true; 571 if (NOT_AN_ID!=hopp) isovalue[hopp] = true; 572 } 573 574 if (!cut && NOT_AN_ID!=fopp && !singular[fopp]) { 575 indir[indir_root(fopp, indir)] = indir_root(f, indir); 576 } 577 } 578 579 FOR(f, m->facets.nb()) { 580 chart[f] = indir_root(f, indir); 581 } 582 } 583 584 585 { // determine quad corner vertices: fill quadcorners attribute 586 FOR (h, m->facet_corners.nb()) { 587 int cnt = 0; 588 index_t cir = h; 589 do { 590 if (NOT_AN_ID==cir) break; 591 cnt += int(isovalue[cir]); 592 cir = fec.next_around_vertex(cir); 593 } while (cir != h); 594 595 if ((NOT_AN_ID==cir && cnt>=2) || cnt>2) { 596 quadcorners[fec.org(h)] = true; 597 } 598 } 599 } 600 601 602 { // fill quadelement attribute 603 604 // this code verifies only if chart boundaries are marked as isovalues + if it has 4 quadcorners 605 // normally it is also necessary to check if the "quad" is a 2d disk (one boundary + Euler characterisic (imagine a torus with a hole)) 606 607 vector<bool> seen(m->facets.nb(), false); 608 FOR(fseed, m->facets.nb()) { 609 if (fseed != chart[fseed]) continue; 610 int nb_quad_corners = 0; 611 std::deque<index_t> Q; 612 Q.push_back(fseed); 613 seen[fseed] = true; 614 bool iso_only_at_boundaries = true; 615 std::vector<index_t> C; 616 while (Q.size()) { 617 index_t f = Q.front(); 618 if (singular[f]) { 619 iso_only_at_boundaries = false; 620 break; 621 } 622 C.push_back(f); 623 Q.pop_front(); 624 FOR(fc, m->facets.nb_corners(f)) { 625 index_t h = m->facets.corner(f, fc); 626 index_t hopp = fec.opposite(h); 627 index_t fopp = NOT_AN_ID==hopp ? NOT_AN_ID : fec.facet(hopp); 628 if (fopp==NOT_AN_ID) { 629 if (isovalue[h] && quadcorners[fec.org(h)]) nb_quad_corners++; 630 } else { 631 if (chart[fopp] != chart[fseed]) { 632 if (isovalue[h]) { 633 if (quadcorners[fec.org(h)]) 634 nb_quad_corners++; 635 continue; 636 } 637 Q.clear(); 638 iso_only_at_boundaries = false; 639 break; 640 } 641 } 642 if (NOT_AN_ID!=fopp && !seen[fopp]) { 643 Q.push_back(fopp); 644 seen[fopp] = true; 645 } 646 } 647 } 648 if (iso_only_at_boundaries && nb_quad_corners == 4) { 649 FOR(i, C.size()) { 650 quadelement[C[i]] = true; 651 } 652 } 653 } 654 } 655 656 { // charts = orig triangles in non-quad regions 657 FOR(f, m->facets.nb()) { 658 if (!quadelement[f]) indir[f] = f; 659 } 660 661 Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid"); 662 FOR(h, m->facet_corners.nb()) { 663 index_t hopp = fec.opposite(h); 664 index_t f = fec.facet(h); 665 index_t fopp = NOT_AN_ID==hopp ? NOT_AN_ID : fec.facet(hopp); 666 if (NOT_AN_ID == fopp || quadelement[f] || quadelement[fopp] || orig_tri_fid[f]!=orig_tri_fid[fopp]) continue; 667 indir[indir_root(fopp, indir)] = indir_root(f, indir); 668 } 669 670 FOR(f, m->facets.nb()) { 671 chart[f] = indir_root(f, indir); 672 } 673 } 674 675 vector<vector<index_t> > v2f = generate_v2f(m); 676 { // mark the boundary between triangles and quads 677 FOR(v, m->vertices.nb()) { 678 TrFan fan = TrFan(v, m, v2f, chart); 679 bool touches_a_quad = false; 680 FOR(i, fan.nb_fan_triangles()) { 681 touches_a_quad = touches_a_quad || quadelement[fan[i].f]; 682 } 683 if (quadcorners[v] || !touches_a_quad || fan.ncharts()<=2) continue; 684 685 index_t first_triangle = NOT_AN_ID; 686 FOR(i, fan.nb_fan_triangles()) { 687 if (quadelement[fan[i].f]) continue; 688 if (NOT_AN_ID==first_triangle) { 689 first_triangle = fan[i].f; 690 } 691 indir[indir_root(first_triangle, indir)] = indir_root(fan[i].f, indir); 692 } 693 } 694 FOR(f, m->facets.nb()) { 695 chart[f] = indir_root(f, indir); 696 } 697 } 698 699 { // fill verts to remove attribute 700 Attribute<bool> verts_to_remove(m->vertices.attributes(), "verts_to_remove"); 701 FOR(v, m->vertices.nb()) { 702 TrFan fan = TrFan(v, m, v2f, chart); 703 bool touches_a_quad = false; 704 FOR(i, fan.nb_fan_triangles()) { 705 touches_a_quad = touches_a_quad || quadelement[fan[i].f]; 706 } 707 if (quadcorners[v]) continue; 708 if ( fan.incomplete_ && fan.ncharts()!=1) continue; 709 if (!fan.incomplete_ && fan.ncharts() >2) continue; 710 if (!touches_a_quad && !added_vertices[v]) continue; 711 verts_to_remove[v] = true; 712 } 713 } 714 } 715 716 /****************************************************************************************************/ try_export_quadtri_from_charts(Mesh * m,vector<BBox> & locked_regions)717 static bool try_export_quadtri_from_charts(Mesh* m, vector<BBox>& locked_regions) { 718 bool modified = false; 719 Attribute<index_t> chart(m->facets.attributes(), "chart"); // TODO document this 720 Attribute<bool> quadelement(m->facets.attributes(), "quadelement"); 721 722 vector<index_t> to_kill(m->facets.nb(), false); 723 724 index_t max_chart_no = 0; 725 FOR(f, m->facets.nb()) { 726 max_chart_no = std::max(max_chart_no, chart[f] + 1); 727 } 728 vector<int> nbtri_in_chart(max_chart_no, 0); 729 FOR(f, m->facets.nb()) { 730 nbtri_in_chart[chart[f]]++; 731 } 732 733 FacetsExtraConnectivity fec(m); 734 index_t nbf = m->facets.nb(); 735 FOR(f, nbf) { 736 geo_assert(3 == m->facets.nb_corners(f)); 737 738 if (nbtri_in_chart[chart[f]] != 2 || !quadelement[f]) continue; 739 FOR(ih, 3) { 740 index_t h = m->facets.corner(f, ih); 741 index_t hopp = fec.opposite(h); 742 if (NOT_AN_ID==hopp) continue; 743 index_t fopp = fec.facet(hopp); 744 geo_assert(NOT_AN_ID != fopp); 745 if (chart[fopp] != chart[f]) continue; 746 if (f < fopp) break; 747 748 vector<index_t> pts; 749 pts.push_back(fec.dest(h)); 750 pts.push_back(fec.dest(fec.next(h))); 751 pts.push_back(fec.org(h)); 752 pts.push_back(fec.dest(fec.next(hopp))); 753 754 bool intersect_locked_region = false; 755 BBox b; 756 FOR(v, 4) { 757 b.add(X(m)[pts[v]]); 758 } 759 FOR(i, locked_regions.size()) { 760 intersect_locked_region = intersect_locked_region || locked_regions[i].intersect(b); 761 } 762 763 if (!intersect_locked_region) { 764 index_t nf = m->facets.create_polygon(pts); 765 modified = true; 766 m->facets.attributes().copy_item(nf ,f); 767 to_kill.push_back(false); 768 to_kill[f] = true; 769 to_kill[fopp] = true; 770 } 771 } 772 } 773 m->facets.delete_elements(to_kill, false); 774 return modified; 775 } 776 777 778 779 sample_triangle(vec3 * ABC,double eps,vector<vec3> & samples)780 static void sample_triangle(vec3 *ABC, double eps, vector<vec3>& samples) { 781 double max_edge_length = 0; 782 FOR(p, 3) max_equal(max_edge_length, (ABC[(p + 1) % 3] - ABC[p]).length()); 783 index_t nb_steps = 10; 784 if (eps>0) min_equal(nb_steps, index_t(max_edge_length / eps + 2)); 785 786 FOR(i, nb_steps)FOR(j, nb_steps - i) { 787 double u = double(i) / double(nb_steps - 1); 788 double v = double(j) / double(nb_steps - 1); 789 samples.push_back(ABC[0] + u*(ABC[1] - ABC[0]) + v*(ABC[2] - ABC[0])); 790 } 791 } 792 //double upper_bound_min_dist2_to_triangles(vec3 P, vector<vec3>& triangles) { 793 // vec3 closest_point; 794 // double l0, l1, l2; 795 // double min_dist2 = 1e20; 796 // FOR(t, triangles.size() / 3) 797 // min_equal(min_dist2, 798 // Geom::point_triangle_squared_distance<vec3>(P, 799 // triangles[t * 3], triangles[t * 3 + 1], triangles[t * 3 + 2], closest_point, l0, l1, l2) 800 // ); 801 // return min_dist2; 802 //} 803 facets_having_a_point_further_than_eps(Mesh * m,Mesh * ref,double epsilon)804 static vector<index_t> facets_having_a_point_further_than_eps(Mesh* m,Mesh* ref,double epsilon) { 805 806 vector<index_t> res; 807 808 vector<BBox> inboxes = facets_bbox(ref); 809 DynamicHBoxes hb; hb.init(inboxes); 810 811 FOR(f, m->facets.nb()) { 812 bool fail = false; 813 vector<vec3> samples; 814 vec3 ABC[3]; 815 FOR(lv,3) ABC[lv] = X(m)[m->facets.vertex(f,lv)]; 816 sample_triangle(ABC, epsilon / 2., samples); 817 if (m->facets.nb_vertices(f) == 4) { 818 FOR(lv, 3) ABC[lv] = X(m)[m->facets.vertex(f, (lv+2)%3)]; 819 sample_triangle(ABC, epsilon / 2., samples); 820 } 821 822 823 FOR(p, samples.size()) { 824 vec3 P = samples[p]; 825 826 double min_dist2 = 1e20; 827 828 BBox bbox; bbox.add(P); bbox.dilate(epsilon/2.); 829 vector<index_t> prim; 830 hb.intersect(bbox, prim); 831 FOR(fid, prim.size()) { 832 index_t other_f = prim[fid]; 833 vec3 closest_point; 834 double l0, l1, l2; 835 min_equal(min_dist2, Geom::point_triangle_squared_distance<vec3>(P, 836 X(ref)[ref->facets.vertex(other_f, 0)], 837 X(ref)[ref->facets.vertex(other_f, 1)], 838 X(ref)[ref->facets.vertex(other_f, 2)], 839 closest_point, l0, l1, l2)); 840 841 if (ref->facets.nb_vertices(other_f) == 4) { 842 min_equal(min_dist2, Geom::point_triangle_squared_distance<vec3>(P, 843 X(ref)[ref->facets.vertex(other_f, 0)], 844 X(ref)[ref->facets.vertex(other_f, 2)], 845 X(ref)[ref->facets.vertex(other_f, 3)], 846 closest_point, l0, l1, l2)); 847 } 848 } 849 850 if (::sqrt(min_dist2) > epsilon / 2.) { 851 fail = true; 852 break; 853 } 854 } 855 if (fail) res.push_back(f); 856 } 857 return res; 858 } 859 860 // Attention, sub-functions of this function need to access to attributes "chart" and "singular" simplify_quad_charts(Mesh * m)861 void simplify_quad_charts(Mesh* m) { 862 std::string msg; 863 if (!surface_is_manifold(m, msg)) plop(msg); 864 865 Mesh m_bak; 866 m_bak.copy(*m); 867 double epsilon = 0;// .4*get_facet_average_edge_size(&m_bak); 868 //epsilon = 0; 869 vector<BBox> locked_regions; 870 // int cnt = 0; 871 for(;;) { 872 vector<index_t> invalid_m ; 873 vector<index_t> invalid_bak ; 874 vector<index_t> intersections; 875 { 876 Attribute<index_t> chart(m->facets.attributes(), "chart"); 877 Attribute<index_t> undo(m->facets.attributes(), "undo"); 878 FOR(fid, m->facets.nb()) { 879 undo[fid] = NOT_AN_ID; 880 } 881 Attribute<bool> verts_to_remove(m->vertices.attributes(), "verts_to_remove"); 882 plop("try_simplify(m, chart, verts_to_remove, undo)"); 883 try_simplify(m, chart, verts_to_remove, undo); 884 plop("try_export_quadtri_from_charts(m, locked_regions)"); 885 try_export_quadtri_from_charts(m, locked_regions); 886 887 plop("check for intersections"); 888 plop(m->facets.nb()); 889 find_self_intersections(m, intersections); 890 FOR(f, intersections.size()) { 891 if (3 == m->facets.nb_vertices(f)) { 892 if (undo[intersections[f]] != NOT_AN_ID) verts_to_remove[undo[intersections[f]]] = false; 893 // GEO::Logger::out("HexDom") << intersections[f] << " " << undo[intersections[f]] << std::endl; 894 } else { 895 geo_assert(4 == m->facets.nb_vertices(f)); 896 BBox inbox; 897 FOR(fv, 4) { 898 inbox.add(X(m)[m->facets.vertex(intersections[f], fv)]); 899 } 900 locked_regions.push_back(inbox); 901 } 902 // Attribute<bool> gna(m->facets.attributes(), "gna"); 903 // gna[intersections[f]] = true; 904 } 905 //intersections.clear();// HACK (simulation de quadhex) 906 907 if (epsilon > 0) { 908 plop("check for Hausdorff distance --- dist to init mesh"); 909 invalid_m = facets_having_a_point_further_than_eps(m, &m_bak, epsilon); 910 plop(invalid_m.size()); 911 FOR(i, invalid_m.size()) { 912 if (3 == m->facets.nb_vertices(invalid_m[i])) { 913 if (undo[invalid_m[i]] != NOT_AN_ID) verts_to_remove[undo[invalid_m[i]]] = false; 914 } 915 else { 916 BBox inbox; 917 FOR(fv, 4) { 918 inbox.add(X(m)[m->facets.vertex(invalid_m[i], fv)]); 919 } 920 locked_regions.push_back(inbox); 921 } 922 } 923 924 plop("check for Hausdorff distance --- dist to new mesh"); 925 invalid_bak = facets_having_a_point_further_than_eps(&m_bak, m, epsilon); 926 plop(invalid_bak.size()); 927 928 FOR(i, invalid_bak.size()) FOR(lv, 3) 929 verts_to_remove[m_bak.facets.vertex(invalid_bak[i], lv)] = false; 930 } 931 } 932 933 934 if (intersections.empty() && invalid_m.empty() && invalid_bak.empty()) break; 935 936 plop("conflict detected"); 937 { 938 Attribute<bool> m_verts_to_remove(m->vertices.attributes(), "verts_to_remove"); 939 Attribute<bool> m_bak_verts_to_remove(m_bak.vertices.attributes(), "verts_to_remove"); 940 geo_assert(m->vertices.nb() == m_bak.vertices.nb()); 941 FOR(v, m->vertices.nb()) m_bak_verts_to_remove[v] = m_verts_to_remove[v]; 942 } 943 944 // char filename[1024], filename2[1024]; 945 // sprintf(filename, "/home/ssloy/tmp/hexdom_nightly/geogram/zdebug%i_bak.geogram", cnt); 946 // sprintf(filename2, "/home/ssloy/tmp/hexdom_nightly/geogram/zdebug%i_simp.geogram", cnt); 947 // cnt++; 948 // mesh_save(m_bak, filename); 949 // mesh_save(*m, filename2); 950 951 m->copy(m_bak); 952 } 953 kill_isolated_vertices(m); 954 } 955 956 957 958 959 } 960 961