1 /* 2 * OGF/Graphite: Geometry and Graphics Programming Library + Utilities 3 * Copyright (C) 2000-2015 INRIA - Project ALICE 4 * 5 * This program is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 2 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program; if not, write to the Free Software 17 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18 * 19 * If you modify this software, you should include a notice giving the 20 * name of the person performing the modification, the date of modification, 21 * and the reason for such modification. 22 * 23 * Contact for Graphite: Bruno Levy - Bruno.Levy@inria.fr 24 * Contact for this Plugin: Nicolas Ray - nicolas.ray@inria.fr 25 * 26 * Project ALICE 27 * LORIA, INRIA Lorraine, 28 * Campus Scientifique, BP 239 29 * 54506 VANDOEUVRE LES NANCY CEDEX 30 * FRANCE 31 * 32 * Note that the GNU General Public License does not permit incorporating 33 * the Software into proprietary programs. 34 * 35 * As an exception to the GPL, Graphite can be linked with the following 36 * (non-GPL) libraries: 37 * Qt, tetgen, SuperLU, WildMagic and CGAL 38 */ 39 40 #include <exploragram/hexdom/PGP.h> 41 #include <exploragram/hexdom/frame.h> 42 #include <exploragram/hexdom/extra_connectivity.h> 43 #include <exploragram/hexdom/geometry.h> 44 #include <exploragram/hexdom/mesh_utils.h> 45 #include <exploragram/hexdom/time_log.h> 46 47 48 49 50 #include <exploragram/hexdom/quadmesher.h> 51 #include <geogram/NL/nl.h> 52 53 #include <algorithm> 54 #include <cmath> 55 #include <queue> 56 57 namespace GEO { 58 PGPopt(Mesh * p_m)59 PGPopt::PGPopt(Mesh* p_m) : m(p_m) { 60 // init fast acces to edges from vertices 61 if (m->edges.nb() == 0) compute_tet_edge_graph(m, v2e, false); 62 else restore_v2e(m, v2e); 63 v2eopp = vector<vector<index_t> >(m->edges.nb(), vector<index_t>()); 64 FOR(e, m->edges.nb()) v2eopp[m->edges.vertex(e, 1)].push_back(e); 65 66 //bind attributes 67 U.bind(m->vertices.attributes(), "U"); 68 B.bind(m->vertices.attributes(), "B"); 69 corr.bind(m->edges.attributes(), "corr"); 70 tij.bind(m->edges.attributes(), "tij"); 71 } 72 is_PGP_singular(index_t c,index_t lf)73 bool PGPopt::is_PGP_singular(index_t c, index_t lf) { 74 index_t v[3]; 75 FOR(i, 3) v[i] = m->cells.facet_vertex(c, lf, i); 76 vec3i t(0, 0, 0); 77 mat3 R; 78 R.load_identity(); 79 FOR(i, 3) { 80 AxisPermutation r = Rij(m, B, v[i], v[(i + 1) % 3]); 81 bool inv; 82 index_t e = edge_from_vertices(v[i], v[(i + 1) % 3], inv); 83 int i1 = int(tij[e][0]); 84 int i2 = int(tij[e][1]); 85 int i3 = int(tij[e][2]); 86 vec3i ltij(i1,i2,i3); 87 t += R*(inv ? -(r.inverse()*ltij) : ltij); 88 R = R* r.inverse().get_mat(); 89 } 90 91 geo_assert(R.is_identity()); 92 return t[0] || t[1] || t[2]; 93 } 94 95 tet_is_PGP_singular_fct(index_t t)96 bool PGPopt::tet_is_PGP_singular_fct(index_t t) { 97 bool is_sing = false; 98 FOR(f, 4) is_sing = is_sing || is_PGP_singular(t, f); 99 return is_sing; 100 } 101 102 face_is_resp(index_t c,index_t lf)103 bool PGPopt::face_is_resp(index_t c, index_t lf) { 104 if (m->cells.adjacent(c, lf) == NO_CELL) return true; 105 return m->cells.adjacent(c, lf) < c; 106 } 107 108 struct TriangleEdgesInSameBasis { try_new_triangleGEO::TriangleEdgesInSameBasis109 bool try_new_triangle(PGPopt* pgp, index_t c, index_t lf) { 110 if (!pgp->face_is_resp(c, lf)) { return false; } // avoid doing twice the same work 111 if (triangle_is_frame_singular(pgp->m, pgp->B, c, lf)) { return false; } // curl correction on singular face is meaningless 112 // without chain basis change 113 index_t vid[3]; // vertices index 114 // init vertex indices 115 FOR(e, 3) vid[e] = pgp->m->cells.facet_vertex(c, lf, e); 116 // init edges 117 FOR(e, 3) { 118 edge[e] = pgp->edge_from_vertices(vid[e], vid[next_mod(e, 3)], inv[e]); 119 geo_assert(edge[e] != NOT_AN_ID); 120 } 121 // express corr of all edges in a common basis with edge_ap[e] 122 FOR(e, 3) { 123 if (inv[e]) edge_ap[e] = Rij(pgp->m, pgp->B, vid[0], vid[next_mod(e, 3)]); 124 else edge_ap[e] = Rij(pgp->m, pgp->B, vid[0], vid[e]); 125 } 126 return true; 127 } 128 129 130 index_t edge[3]; 131 bool inv[3]; 132 AxisPermutation edge_ap[3]; 133 }; 134 135 136 137 138 // ___ __ 139 // / _ \ _ _ ___ / _|___ _ _ _ __ ___ 140 // | (_) | ' \/ -_) | _/ _ \ '_| ' \(_-< 141 // \___/|_||_\___| |_| \___/_| |_|_|_/__/ 142 // 143 144 145 wish_angle_corr(index_t e,bool inv)146 vec3 PGPopt::wish_angle_corr(index_t e, bool inv) { 147 vec3 c(0, 0, 0); 148 if (corr.is_bound()) { // CubeCover ne l'utilise pas forcement 149 c = corr[e]; 150 if (inv) { 151 AxisPermutation chg = Rij(m, B, m->edges.vertex(e, 0), m->edges.vertex(e, 1)); 152 c = -(chg.inverse() *c); 153 } 154 } 155 return c; 156 } 157 wish_angle_edge_geom(index_t e,bool inv)158 vec3 PGPopt::wish_angle_edge_geom(index_t e, bool inv) { 159 index_t org = m->edges.vertex(e, 0); 160 index_t dest = m->edges.vertex(e, 1); 161 if (inv) std::swap(org, dest); 162 163 AxisPermutation ap = Rij(m, B, org, dest); 164 165 mat3 frame = B[org] + Frame(B[dest]).apply_permutation(ap); 166 frame *= 0.5; 167 frame = invert_columns_norm(frame); 168 vec3 angle = frame.transpose() * (m->vertices.point(dest) - m->vertices.point(org)); 169 return 2.*M_PI *angle; // one cycle length is edgelength_ 170 } 171 wish_angle(index_t e,bool inv)172 vec3 PGPopt::wish_angle(index_t e, bool inv) { 173 return wish_angle_edge_geom(e, inv) + PGPopt::wish_angle_corr(e, inv); 174 } 175 176 177 178 179 180 // _ _ _ 181 // ___ _ __| |_(_)_ __ (_)______ __ ___ _ _ _ _ 182 // / _ \ '_ \ _| | ' \| |_ / -_) / _/ _ \ '_| '_| 183 // \___/ .__/\__|_|_|_|_|_/__\___|_\__\___/_| |_| 184 // |_| 185 186 187 /////////////////////////////////////////////////// 188 //optimize_corr 189 /////////////////////////////////////////////////// 190 // problem: the objective one form derived from the frame field is not close, leading to PGP singularities (T_junctions) 191 // solution: this function computes a one form such that : 192 // *its norm is minimal 193 // *adding it to the original objective one form, make it close... everywhere but on frame field singularities 194 // param : max_corr_prop typically in [0,.5] (0 => no correction, .5 => take initial field into account for only 50% 195 // 196 // note: it is the same energy as in cubecover, exept that there is no integer constraints, and the frame singularities are ignored 197 198 optimize_corr(double max_corr_prop)199 void PGPopt::optimize_corr(double max_corr_prop) { 200 FOR(e, m->edges.nb()) corr[e] = vec3(0, 0, 0); 201 if (max_corr_prop == 0) { 202 return; 203 } 204 cubcover(true); 205 // clamp the result: very naive way to avoid too large corrections. May be improved. 206 FOR(e, m->edges.nb()) { 207 // BasisChg chg = edge_basis_change(e, false); 208 vec3 geom_form = wish_angle_edge_geom(e, false); 209 vec3 corr_form = wish_angle_corr(e, false); 210 double scale = 1.; 211 double cl = corr_form.length(); 212 double gl = geom_form.length(); 213 if (cl > max_corr_prop* gl) scale = max_corr_prop* gl / cl; 214 FOR(d, 3) corr[e][d] = scale * corr[e][d]; 215 } 216 } 217 218 219 move_U_to_corner()220 void PGPopt::move_U_to_corner() { 221 Attribute<vec3> UC(m->cell_corners.attributes(), "U"); 222 Attribute<bool> has_param(m->cell_facets.attributes(), "has_param"); 223 224 FOR(c, m->cells.nb()) { 225 226 bool all_faces_have_param = true; 227 // flag has param 228 FOR(lf, 4) { 229 index_t f = m->cells.facet(c, lf); 230 has_param[f] = true; 231 if (triangle_is_frame_singular(m, B, c, lf))has_param[f] = false; 232 else if (is_PGP_singular(c, lf)) has_param[f] = false; 233 all_faces_have_param = all_faces_have_param && has_param[f]; 234 } 235 236 237 index_t org = 0; 238 239 // find a seed that will properly reconstruct all valid triangle param 240 if (!all_faces_have_param) { 241 bool can_be_ref[4] = { true, true, true, true }; 242 FOR(lf, 4) { 243 index_t f = m->cells.facet(c, lf); 244 if (has_param[f]) { 245 bool local_can_be_ref[4] = { false, false, false, false }; 246 FOR(lv, 3) local_can_be_ref[m->cells.descriptor(c).facet_vertex[lf][lv]] = true; 247 FOR(i, 4) can_be_ref[i] = can_be_ref[i] && local_can_be_ref[i]; 248 } 249 } 250 org = NOT_AN_ID; 251 FOR(i, 4) if (can_be_ref[i]) org = i; 252 } 253 254 255 if (org != NOT_AN_ID) FOR(i, 4) { 256 index_t corner = m->cells.corner(c, i); 257 UC[corner] = U[m->cells.vertex(c, i)]; 258 if (i != org) { 259 AxisPermutation change = Rij(m, B, m->cells.vertex(c, org), m->cells.vertex(c, i)); 260 UC[corner] = change.inverse() * U[m->cells.vertex(c, i)]; 261 bool inv; 262 index_t e = edge_from_vertices(m->cells.vertex(c, org), m->cells.vertex(c, i), inv); 263 geo_assert(e != NOT_AN_ID); 264 vec3 t2 = !inv ? tij[e] : -(change.inverse()*tij[e]); 265 UC[corner] += vec3(t2[0], t2[1], t2[2]); 266 267 } 268 } 269 else geo_assert(!has_param[m->cells.facet(c, 0)] && !has_param[m->cells.facet(c, 1)] && !has_param[m->cells.facet(c, 2)] && !has_param[m->cells.facet(c, 3)]); 270 } 271 272 // TODO REMOVE THAT !!!!! just to TEST 273 //FOR(c, m->cell_corners.nb()) UC[c] = 2*UC[c]; 274 } 275 276 /* 277 * _ _ _ ___ ___ ___ 278 * ___ _ __| |_(_)_ __ (_)______ | _ \/ __| _ \ 279 * / _ \ '_ \ _| | ' \| |_ / -_) | _/ (_ | _/ 280 * \___/ .__/\__|_|_|_|_|_/__\___|_|_| \___|_| 281 * |_| 282 */ 283 optimize_PGP()284 void PGPopt::optimize_PGP() { 285 Attribute<vec3> lockU(m->vertices.attributes(), "lockU"); 286 // Create and initialize OpenNL context 287 nlNewContext(); 288 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); 289 nlSolverParameteri(NL_NB_VARIABLES, NLint(6 * m->vertices.nb())); 290 291 nlBegin(NL_SYSTEM); 292 FOR(v, m->vertices.nb())FOR(d, 3) if (std::fabs(lockU[v][d]) > 0)FOR(c, 2) { 293 nlSetVariable(6 * v + 2 * d + c, 1 - c); 294 nlLockVariable(6 * v + 2 * d + c); 295 } 296 297 nlBegin(NL_MATRIX); 298 FOR(e, m->edges.nb()) { 299 AxisPermutation ap = Rij(m, B, m->edges.vertex(e, 0), m->edges.vertex(e, 1)); 300 301 vec3 theta = wish_angle(e, false); 302 FOR(d, 3) { 303 double c = cos(theta[d]); 304 double s = sin(theta[d]); 305 index_t off0 = 6 * m->edges.vertex(e, 0) + 2 * d; 306 307 nlBegin(NL_ROW); 308 FOR(dd, 3) if (ap.get_mat()(dd, d) != 0) 309 nlCoefficient(6 * m->edges.vertex(e, 1) + 2 * dd, -1.); 310 nlCoefficient(off0, c); 311 nlCoefficient(off0 + 1, s); 312 nlEnd(NL_ROW); 313 nlBegin(NL_ROW); 314 FOR(dd, 3) 315 nlCoefficient(6 * m->edges.vertex(e, 1) + 2 * dd + 1, -ap.get_mat()(dd, d)); 316 nlCoefficient(off0, -s); 317 nlCoefficient(off0 + 1, c); 318 nlEnd(NL_ROW); 319 } 320 } 321 322 nlEnd(NL_MATRIX); 323 nlEnd(NL_SYSTEM); 324 // Solve and get solution 325 nlSolve(); 326 327 FOR(v, m->vertices.nb()) FOR(d, 3) 328 U[v][d] = (.5 / M_PI) * atan2(nlGetVariable(6 * v + 2 * d + 1), nlGetVariable(6 * v + 2 * d)); 329 330 nlDeleteContext(nlGetCurrent()); 331 332 333 FOR(e, m->edges.nb()) { 334 index_t i = m->edges.vertex(e, 0); 335 index_t j = m->edges.vertex(e, 1); 336 AxisPermutation rij = Rij(m, B, i, j); 337 338 vec3 gij = wish_angle(e, false) / (2.*M_PI); 339 340 FOR(d, 3) { 341 tij[e][d] = int(round(-(rij.inverse().get_mat()*U[j])[d] - gij[d] + U[i][d])); 342 } 343 344 } 345 move_U_to_corner(); 346 } 347 348 349 350 351 352 353 // _____ _ _____ 354 // / ____| | | / ____| 355 // | | _ _| |__ ___| | _____ _____ _ __ 356 // | | | | | | '_ \ / _ \ | / _ \ \ / / _ \ '__| 357 // | |___| |_| | |_) | __/ |___| (_) \ V / __/ | 358 // \_____\__,_|_.__/ \___|\_____\___/ \_/ \___|_| 359 // 360 // 361 struct TetHalfedge { TetHalfedgeGEO::TetHalfedge362 TetHalfedge(index_t p_cell, index_t p_org, index_t p_dest) { 363 cell = p_cell; org = p_org; dest = p_dest; 364 } 365 index_t cell; 366 index_t org; 367 index_t dest; 368 }; get_non_nulledge(index_t c,index_t cf,Attribute<bool> & nulledge)369 index_t PGPopt::get_non_nulledge(index_t c, index_t cf, Attribute<bool>& nulledge){ 370 int nbzero = 0; 371 FOR(cfv, 3) { 372 bool inv; 373 index_t e = edge_from_vertices(m->cells.facet_vertex(c, cf, cfv), m->cells.facet_vertex(c, cf, next_mod(cfv, 3)), inv); 374 geo_assert(e != NOT_AN_ID); 375 if (nulledge[e]) nbzero++; 376 } 377 if (nbzero == 2 && !triangle_is_frame_singular(m, B, c, cf)) { 378 FOR(cfv, 3) { 379 bool inv; 380 index_t e = edge_from_vertices(m->cells.facet_vertex(c, cf, cfv), m->cells.facet_vertex(c, cf, next_mod(cfv, 3)), inv); 381 geo_assert(e != NOT_AN_ID); 382 if (!nulledge[e]) { 383 nulledge[e] = true; 384 return cfv; 385 } 386 } 387 } 388 return NOT_AN_ID; 389 } 390 391 392 393 mark_null_edges(Attribute<bool> & nulledge)394 void PGPopt::mark_null_edges(Attribute<bool>& nulledge) { 395 396 Attribute<vec3> lockU(m->vertices.attributes(), "lockU"); 397 index_t alot = std::numeric_limits<index_t>::max() - 10; 398 FOR(e, m->edges.nb()) nulledge[e] = false; 399 400 vector<index_t> offset_from_org; 401 vector<index_t> dest; 402 403 // compute covering tree 404 cell_edges_in_RCS(m, offset_from_org, dest); 405 std::vector<index_t > dist(m->vertices.nb(), alot); 406 407 { 408 FOR(seed, m->vertices.nb()) { 409 if (dist[seed] != alot) continue; 410 dist[seed] = 0; 411 std::queue<index_t> queue; 412 queue.push(seed); 413 while (!queue.empty()) { 414 index_t v = queue.front(); 415 queue.pop(); 416 for (index_t ed = offset_from_org[v]; ed < offset_from_org[v + 1]; ed++) { 417 index_t vopp = dest[ed]; 418 if (dist[vopp] == alot) { 419 bool inv; 420 index_t e = edge_from_vertices(v, vopp, inv); 421 geo_assert(e != NOT_AN_ID); 422 nulledge[e] = true; 423 dist[vopp] = dist[v] + 1; 424 queue.push(vopp); 425 } 426 } 427 } 428 } 429 } 430 431 std::queue<TetHalfedge> queue; 432 // init 433 FOR(c, m->cells.nb()) FOR(cf, 4) { 434 index_t cfe = get_non_nulledge(c, cf, nulledge); 435 if (cfe != NOT_AN_ID) 436 queue.push(TetHalfedge(c, m->cells.facet_vertex(c, cf, cfe), m->cells.facet_vertex(c, cf, next_mod(cfe, 3)))); 437 } 438 while (!queue.empty()) { 439 TetHalfedge cur = queue.front(); 440 queue.pop(); 441 index_t cir = cur.cell; 442 do { 443 FOR(cf, 4) { 444 index_t cfe = get_non_nulledge(cir, cf, nulledge); 445 if (cfe != NOT_AN_ID) 446 queue.push(TetHalfedge(cir, m->cells.facet_vertex(cir, cf, cfe), m->cells.facet_vertex(cir, cf, next_mod(cfe, 3)))); 447 } 448 cir = next_cell_around_oriented_edge(m, cir, cur.org, cur.dest); 449 450 if (cir == NOT_AN_ID) { 451 cir = cur.cell; 452 while (next_cell_around_oriented_edge(m, cir, cur.dest, cur.org) != NOT_AN_ID) 453 cir = next_cell_around_oriented_edge(m, cir, cur.dest, cur.org); 454 } 455 } while (cir != cur.cell); 456 } 457 458 459 } 460 461 // // _____ _ _____ 462 // // / ____| | | / ____| 463 // // | | _ _| |__ ___| | _____ _____ _ __ 464 // // | | | | | | '_ \ / _ \ | / _ \ \ / / _ \ '__| 465 // // | |___| |_| | |_) | __/ |___| (_) \ V / __/ | 466 // // \_____\__,_|_.__/ \___|\_____\___/ \_/ \___|_| 467 // // 468 // // 469 grow_ball(Attribute<bool> & tet_in_ball)470 void PGPopt::grow_ball(Attribute<bool>& tet_in_ball) { 471 472 vector<bool> touched(m->vertices.nb(), false); 473 vector<bool> tet_is_singular(m->cells.nb(), false); 474 FOR(c, m->cells.nb()) FOR(lf, 4) tet_is_singular[c] = tet_is_singular[c] || triangle_is_frame_singular(m, B, c, lf); 475 476 FOR(c, m->cells.nb()) tet_in_ball[c] = false; 477 FOR(seed, m->cells.nb()) { 478 if (tet_in_ball[seed]) continue; 479 if (tet_is_singular[seed]) continue; 480 if (touched[m->cells.vertex(seed, 0)]) continue; 481 if (touched[m->cells.vertex(seed, 1)]) continue; 482 if (touched[m->cells.vertex(seed, 2)]) continue; 483 if (touched[m->cells.vertex(seed, 3)]) continue; 484 485 tet_in_ball[seed] = true; 486 FOR(lv, 4) touched[m->cells.vertex(seed, lv)] = true; 487 488 std::queue<index_t> queue; 489 queue.push(seed); 490 while (!queue.empty()) { 491 index_t c = queue.front(); 492 queue.pop(); 493 FOR(lf, 4) { 494 index_t copp = m->cells.adjacent(c, lf); 495 if (copp == NOT_AN_ID) continue; 496 if (tet_in_ball[copp]) continue; 497 if (tet_is_singular[copp]) continue; 498 499 // try to reach a new vertex 500 bool all_are_fine = true; 501 FOR(lvopp, 4) { 502 index_t vopp = m->cells.vertex(copp, lvopp); 503 bool is_in_c = false; 504 FOR(lv, 4) is_in_c = is_in_c || (vopp == m->cells.vertex(c, lv)); 505 all_are_fine = all_are_fine && (is_in_c || !touched[vopp]); 506 } 507 // try to glue to two faces 508 int nb_done_neigs = 0; 509 FOR(lfopp, 4) 510 if (m->cells.adjacent(copp, lfopp) != NOT_AN_ID) 511 if (tet_in_ball[m->cells.adjacent(copp, lfopp)]) 512 nb_done_neigs++; 513 if (!all_are_fine && nb_done_neigs ==1) continue; 514 515 FOR(lvopp, 4) touched[m->cells.vertex(copp, lvopp)] = true; 516 tet_in_ball[copp] = true; 517 queue.push(copp); 518 } 519 } 520 } 521 Attribute<bool> null_edge(m->edges.attributes(), "null"); 522 FOR(e, m->edges.nb()) null_edge[e] = false; 523 FOR(c, m->cells.nb()) if (tet_in_ball[c]) { 524 bool inv; 525 FOR(lv, 3)for (index_t lv2 = lv + 1; lv2 < 4; lv2++) FOR(d, 3) 526 null_edge[edge_from_vertices(m->cells.vertex(c, lv), m->cells.vertex(c, lv2), inv)] = true; 527 } 528 } 529 530 cubcover(bool compute_only_corr)531 void PGPopt::cubcover(bool compute_only_corr) { 532 geo_argused(compute_only_corr); 533 Logger::warn("PGP") << "Cubecover not implemented in the public version" 534 << std::endl; 535 } 536 537 } 538