| 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 <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 | |||
| 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 | |||
| 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 | |||
| 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 | |||
| 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 { | ||
| 109 | 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 | |||
| 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 | |||
| 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 | |||
| 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 | |||
| 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 | |||
| 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 | |||
| 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 { | ||
| 362 | ✗ | 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 | }; | ||
| 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 | |||
| 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 | |||
| 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 | |||
| 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 |