| 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/polygon.h> | ||
| 41 | #include <geogram/NL/nl.h> | ||
| 42 | namespace GEO { | ||
| 43 | |||
| 44 | ✗ | vec2 Poly2d::barycenter() { | |
| 45 | ✗ | vec2 bary(0, 0); | |
| 46 | ✗ | FOR(fv, index_t(pts.size())) { | |
| 47 | ✗ | bary = bary + (1. / double(pts.size()))*pts[fv]; | |
| 48 | } | ||
| 49 | ✗ | return bary; | |
| 50 | } | ||
| 51 | |||
| 52 | static int dump_contour_save_id = 0; | ||
| 53 | ✗ | void Poly2d::dump_contour() { | |
| 54 | ✗ | index_t nbv = pts.size(); | |
| 55 | ✗ | Mesh export_mesh; | |
| 56 | ✗ | export_mesh.vertices.create_vertices(nbv); | |
| 57 | ✗ | FOR(i, nbv) X(&export_mesh)[i] = vec3(pts[i][0], pts[i][1], 0); | |
| 58 | ✗ | vector<index_t> num; | |
| 59 | ✗ | FOR(i, nbv) num.push_back(i); | |
| 60 | ✗ | export_mesh.facets.create_polygon(num); | |
| 61 | char filename[1024]; | ||
| 62 | ✗ | snprintf(filename, 1024, "C:/DATA/2dcontours/nimp2D%i.obj", dump_contour_save_id++); | |
| 63 | ✗ | mesh_save(export_mesh, filename); | |
| 64 | ✗ | } | |
| 65 | |||
| 66 | ✗ | void Poly3d::dump_contour() { | |
| 67 | ✗ | index_t nbv = pts.size(); | |
| 68 | ✗ | Mesh export_mesh; | |
| 69 | ✗ | export_mesh.vertices.create_vertices(nbv); | |
| 70 | ✗ | FOR(i, nbv) X(&export_mesh)[i] = vec3(pts[i][0], pts[i][1], pts[i][2] ); | |
| 71 | ✗ | vector<index_t> num; | |
| 72 | ✗ | FOR(i, nbv) num.push_back(i); | |
| 73 | ✗ | export_mesh.facets.create_polygon(num); | |
| 74 | char filename[1024]; | ||
| 75 | ✗ | snprintf(filename, 1024, "C:/DATA/2dcontours/nimp3D%i.obj", dump_contour_save_id++); | |
| 76 | ✗ | mesh_save(export_mesh, filename); | |
| 77 | ✗ | } | |
| 78 | |||
| 79 | // returns 1024. if concave angle is encountered or if proposed triangle contains one of pts | ||
| 80 | // otherwise returns max angle of the proposed triangle | ||
| 81 | ✗ | double Poly2d::cost(index_t i, index_t j, index_t k) { | |
| 82 | ✗ | vec2 C[3] = { pts[i], pts[j], pts[k] }; | |
| 83 | ✗ | double m = 0; | |
| 84 | ✗ | FOR(v, 3) { | |
| 85 | // note that angle is not the angle inside the triangle, but its complement | ||
| 86 | // angle variable has the "direction" information, thus it is negative for concave angles (right turn) and positive for convex angles (left turn) | ||
| 87 | ✗ | double angle = atan2( | |
| 88 | ✗ | det(C[(v + 1) % 3] - C[(v + 0) % 3], C[(v + 2) % 3] - C[(v + 1) % 3]), | |
| 89 | ✗ | dot(C[(v + 1) % 3] - C[(v + 0) % 3], C[(v + 2) % 3] - C[(v + 1) % 3]) | |
| 90 | ); | ||
| 91 | ✗ | if (angle <= 0) return 1024.; | |
| 92 | ✗ | m = std::max(m, M_PI - angle); | |
| 93 | } | ||
| 94 | |||
| 95 | ✗ | FOR(other, pts.size()) { // TODO c'est con de faire ça, vaut mieux regarder si le triangle est inversé (ça ne gere pas tout [comme, d'ailleurs, le teste courant!]) | |
| 96 | ✗ | if (other == i || other == j || other == k) continue; | |
| 97 | ✗ | vec2 P = pts[other]; | |
| 98 | ✗ | bool inside = true; | |
| 99 | ✗ | FOR(l, 3) { | |
| 100 | ✗ | inside = inside && (det(normalize(C[(l + 1) % 3] - C[l]), normalize(P - C[l])) > 0); | |
| 101 | } | ||
| 102 | ✗ | if (inside) return 1024.; | |
| 103 | } | ||
| 104 | ✗ | return m; | |
| 105 | } | ||
| 106 | |||
| 107 | // this function has O(n^4) computational cost | ||
| 108 | ✗ | bool Poly2d::try_triangulate_minweight(vector<index_t>& triangles) { | |
| 109 | ✗ | triangles.clear(); | |
| 110 | ✗ | index_t n = pts.size(); | |
| 111 | ✗ | geo_assert(n >= 3); | |
| 112 | |||
| 113 | //if (n == 3) { | ||
| 114 | // FOR(v, 3) { | ||
| 115 | // triangles.push_back(v); | ||
| 116 | // } | ||
| 117 | // return true; | ||
| 118 | //} | ||
| 119 | |||
| 120 | // we store in this table results of subproblems | ||
| 121 | // table[i*n + j] stores the triangulation cost for points from i to j | ||
| 122 | // the entry table[0*n + n-1] has the final result. | ||
| 123 | ✗ | std::vector<double> table(n*n, 0.); | |
| 124 | |||
| 125 | // this table stores triangle indices: for each subproblem (i,j) we have table[i*n + j]==k, i.e. the triangle is (i,k,j) | ||
| 126 | ✗ | std::vector<index_t> tri(n*n, index_t(-1)); | |
| 127 | |||
| 128 | // note that the table is filled in diagonals; elements below main diagonal are not used at all | ||
| 129 | ✗ | for (index_t pbsize = 2; pbsize < n; pbsize++) { | |
| 130 | ✗ | for (index_t i = 0, j = pbsize; j < n; i++, j++) { | |
| 131 | // recall that we are testing triangle (i,k,j) which splits the problem (i,j) into | ||
| 132 | // two smaller subproblems (i,k) and (k,j) | ||
| 133 | ✗ | double minv = 1e20; | |
| 134 | ✗ | index_t mink = index_t(-1); | |
| 135 | ✗ | for (index_t k = i + 1; k < j; k++) { | |
| 136 | ✗ | double val = table[i*n + k] + table[k*n + j] + cost(i, k, j); | |
| 137 | ✗ | if (minv <= val) continue; | |
| 138 | ✗ | minv = val; | |
| 139 | ✗ | mink = k; | |
| 140 | } | ||
| 141 | ✗ | geo_assert(mink!=index_t(-1)); | |
| 142 | ✗ | table[i*n + j] = minv; | |
| 143 | ✗ | tri[i*n + j] = mink; | |
| 144 | } | ||
| 145 | } | ||
| 146 | |||
| 147 | // if (table[n-1] >= 1024.) return false; | ||
| 148 | |||
| 149 | ✗ | vector<index_t> Q(1, n - 1); | |
| 150 | ✗ | FOR(t, Q.size()) { | |
| 151 | ✗ | index_t idx = Q[t]; | |
| 152 | |||
| 153 | ✗ | index_t i = idx / n; | |
| 154 | ✗ | index_t k = tri[idx]; | |
| 155 | ✗ | index_t j = idx % n; | |
| 156 | |||
| 157 | ✗ | geo_assert(i!=index_t(-1) && k != index_t(-1) && j!=index_t(-1)); | |
| 158 | ✗ | triangles.push_back(i); | |
| 159 | ✗ | triangles.push_back(k); | |
| 160 | ✗ | triangles.push_back(j); | |
| 161 | |||
| 162 | ✗ | if (k + 2 <= j) Q.push_back(k*n + j); | |
| 163 | ✗ | if (i + 2 <= k) Q.push_back(i*n + k); | |
| 164 | } | ||
| 165 | |||
| 166 | ✗ | if (table[n - 1] >= 1024.) { | |
| 167 | ✗ | plop("may dump_contour for debug...");//dump_contour(); | |
| 168 | } | ||
| 169 | ✗ | return table[n-1] < 1024.; | |
| 170 | ✗ | } | |
| 171 | |||
| 172 | |||
| 173 | // find parity of original points | ||
| 174 | ✗ | index_t Poly2d::parity_of_original_points() { | |
| 175 | ✗ | index_t nbv = pts.size(); | |
| 176 | ✗ | index_t dec = 0; | |
| 177 | ✗ | double dec_score[2] = { 0, 0 }; | |
| 178 | ✗ | FOR(q, nbv / 2) { | |
| 179 | ✗ | FOR(d, 2) | |
| 180 | ✗ | dec_score[d] = std::max(dec_score[d], std::abs( | |
| 181 | ✗ | det(pts[(q * 2 + 0 + d) % nbv] - pts[(q * 2 + 1 + d) % nbv], | |
| 182 | ✗ | pts[(q * 2 + 2 + d) % nbv] - pts[(q * 2 + 1 + d) % nbv]))); | |
| 183 | } | ||
| 184 | ✗ | if (dec_score[1] < dec_score[0])dec = 1; | |
| 185 | ✗ | return dec; | |
| 186 | } | ||
| 187 | |||
| 188 | |||
| 189 | ✗ | bool Poly2d::middle_point_quadrangulate(vector<index_t>& quads) { | |
| 190 | ✗ | index_t nbv = pts.size(); | |
| 191 | ✗ | vec2 G = barycenter(); | |
| 192 | ✗ | index_t dec = parity_of_original_points(); | |
| 193 | ✗ | FOR(q, nbv / 2) { | |
| 194 | ✗ | quads.push_back(nbv); | |
| 195 | ✗ | FOR(v, 3) quads.push_back((q * 2 + 1 - dec + v) % nbv); | |
| 196 | } | ||
| 197 | ✗ | pts.push_back(G); | |
| 198 | ✗ | return true; | |
| 199 | } | ||
| 200 | |||
| 201 | |||
| 202 | ✗ | bool Poly2d::quads_are_valid(vector<index_t>& quads) { | |
| 203 | |||
| 204 | // geometric criteria | ||
| 205 | ✗ | FOR(q, quads.size() / 4) { | |
| 206 | ✗ | FOR(e, 4) { | |
| 207 | ✗ | vec2 v0 = normalize(pts[quads[4 * q + next_mod(e, 4)]] - pts[quads[4 * q + e]]); | |
| 208 | ✗ | vec2 v1 = normalize(pts[quads[4 * q + prev_mod(e, 4)]] - pts[quads[4 * q + e]]); | |
| 209 | ✗ | if (det(v0, v1) < sin(M_PI / 8.)) return false; | |
| 210 | } | ||
| 211 | } | ||
| 212 | ✗ | return true; | |
| 213 | } | ||
| 214 | |||
| 215 | |||
| 216 | |||
| 217 | struct Contour2D { | ||
| 218 | ✗ | void resize(index_t n) { pos_.resize(n); angu_.resize(n);vid_.resize(n);} | |
| 219 | ✗ | void compute_angu() { | |
| 220 | //plop("compute_angu() will not sffice to capture sing 5"); | ||
| 221 | ✗ | angu_.resize(pos_.size()); | |
| 222 | ✗ | FOR(v, pos_.size()) { | |
| 223 | ✗ | angu_[v] = 1; | |
| 224 | ✗ | vec2 d0 = pos(v) - pos(int(v)-1); | |
| 225 | ✗ | vec2 d1 = pos(v+1) - pos(v); | |
| 226 | ✗ | double angle = atan2(det(d0, d1), dot(d0, d1)); | |
| 227 | ✗ | if (angle < M_PI / 4.) angu_[v] = 0; | |
| 228 | ✗ | if (angle < -M_PI / 4.) angu_[v] = -1; | |
| 229 | } | ||
| 230 | ✗ | } | |
| 231 | |||
| 232 | vec2 normal(int v) {// equals 0 for the singularity | ||
| 233 | mat2 R90; R90(0, 0) = 0; R90(0, 1) = -1; R90(1, 0) = 1; R90(1, 1) = 0; | ||
| 234 | return normalize(R90*(pos(v+1) - pos(v - 1))); | ||
| 235 | } | ||
| 236 | |||
| 237 | ✗ | vec2& pos(int i) { return aupp(i, pos_); } | |
| 238 | ✗ | int& angu(int i) { return aupp(i, angu_); } | |
| 239 | ✗ | int& vid(int i) { return aupp(i, vid_); } | |
| 240 | |||
| 241 | ✗ | vec2& pos(index_t i) { return aupp(i, pos_); } | |
| 242 | ✗ | int& angu(index_t i) { return aupp(i, angu_); } | |
| 243 | ✗ | int& vid(index_t i) { return aupp(i, vid_); } | |
| 244 | |||
| 245 | |||
| 246 | void show() { | ||
| 247 | GEO::Logger::out("HexDom") << "\npos.size = " << pos_.size() << std::endl; FOR(i, pos_.size()) std::cerr << pos_[i] << "\t"; | ||
| 248 | GEO::Logger::out("HexDom") << "\nangu.size = " << angu_.size() << std::endl; FOR(i, angu_.size()) std::cerr << angu_[i] << "\t"; | ||
| 249 | GEO::Logger::out("HexDom") << "\nvid.size = " << vid_.size() << std::endl; FOR(i, vid_.size()) std::cerr << vid_[i] << "\t"; | ||
| 250 | } | ||
| 251 | |||
| 252 | ✗ | void remove(int i) { | |
| 253 | ✗ | i = i%int(pos_.size()); | |
| 254 | ✗ | pos_.erase(pos_.begin() + i); | |
| 255 | ✗ | angu_.erase(angu_.begin() + i); | |
| 256 | ✗ | vid_.erase(vid_.begin() + i); | |
| 257 | |||
| 258 | ✗ | } | |
| 259 | |||
| 260 | |||
| 261 | vector<vec2> pos_; | ||
| 262 | vector<int> angu_; | ||
| 263 | vector<int> vid_; | ||
| 264 | }; | ||
| 265 | |||
| 266 | static int export_debug_mesh_id = 0; | ||
| 267 | struct QuadrangulateWithOneSingularity { | ||
| 268 | ✗ | QuadrangulateWithOneSingularity(vector<vec2>& p_pts, vector<index_t>& p_quads) | |
| 269 | ✗ | :pts(p_pts), quads(p_quads) { | |
| 270 | ✗ | R90(0, 0) = 0; R90(0, 1) = -1; R90(1, 0) = 1; R90(1, 1) = 0; | |
| 271 | ✗ | } | |
| 272 | |||
| 273 | // returns the index of the singularity | ||
| 274 | ✗ | int init_contour(vector<int>& angu,int sing_valence) { | |
| 275 | ✗ | index_t offset = 0; | |
| 276 | // find the best offset | ||
| 277 | ✗ | double best_dist2 = 1e20; | |
| 278 | ✗ | contour.resize(pts.size()+1); | |
| 279 | ✗ | contour.pos(0) = vec2(0, 0); | |
| 280 | ✗ | FOR(off, pts.size()) { | |
| 281 | ✗ | vec2 dir(1, 0); | |
| 282 | ✗ | FOR(v, pts.size()) { | |
| 283 | ✗ | contour.pos(v + 1) = contour.pos(v) + dir; | |
| 284 | ✗ | if (aupp(off + v + 1, angu) < 0) dir = -(R90*dir); | |
| 285 | ✗ | if (aupp(off + v + 1, angu) > 0) dir = R90*dir; | |
| 286 | } | ||
| 287 | ✗ | vec2 diag = contour.pos(-1) ; | |
| 288 | |||
| 289 | ✗ | if (sing_valence == 3 && diag.x == -diag.y && diag.length2() < best_dist2) { | |
| 290 | ✗ | best_dist2 = diag.length2(); | |
| 291 | ✗ | offset = off; | |
| 292 | }; | ||
| 293 | ✗ | if (sing_valence == 5 && diag.x == diag.y && diag.length2() < best_dist2) { | |
| 294 | ✗ | best_dist2 = diag.length2(); | |
| 295 | ✗ | offset = off; | |
| 296 | }; | ||
| 297 | } | ||
| 298 | // decal gridpos w.r.t offset | ||
| 299 | ✗ | vec2 dir(1, 0); | |
| 300 | ✗ | FOR(v, pts.size()) { | |
| 301 | ✗ | contour.pos(v + 1) = contour.pos(v) + dir; | |
| 302 | ✗ | if (aupp(offset + v + 1, angu) < 0) dir = -(R90*dir); | |
| 303 | ✗ | if (aupp(offset + v + 1, angu) > 0) dir = R90*dir; | |
| 304 | } | ||
| 305 | |||
| 306 | // define mapping contour -> pts | ||
| 307 | ✗ | FOR(v, pts.size() ) contour.vid(v ) = int(offset + v ) % int(pts.size()); | |
| 308 | ✗ | contour.vid_.back() = contour.vid_.front(); | |
| 309 | |||
| 310 | // singularity on border TODO CHECK angu on singularity ! | ||
| 311 | ✗ | if ((contour.pos(-1) - contour.pos(0)).length2() < .1) { contour.remove(int(contour.pos_.size()) - 1); contour.compute_angu(); return 0; } | |
| 312 | |||
| 313 | // add pts | ||
| 314 | ✗ | vec2 A = contour.pos(0); | |
| 315 | ✗ | vec2 B = contour.pos(-1); | |
| 316 | ✗ | if (sing_valence == 3) for (int i = int(B.x + 1.0); i < int(A.x); i++) { | |
| 317 | ✗ | contour.pos_.push_back(vec2(i, B.y)); | |
| 318 | ✗ | contour.vid_.push_back(int(pts.size())); | |
| 319 | ✗ | pts.push_back(contour.pos_.back()); | |
| 320 | } | ||
| 321 | |||
| 322 | ✗ | if (sing_valence == 5) for (int i = int(B.x - 1); i > int(A.x); i--) { | |
| 323 | ✗ | contour.pos_.push_back(vec2(double(i), B.y)); | |
| 324 | ✗ | contour.vid_.push_back(int(pts.size())); | |
| 325 | ✗ | pts.push_back(contour.pos_.back()); | |
| 326 | } | ||
| 327 | |||
| 328 | ✗ | index_t singularity_index = contour.pos_.size(); | |
| 329 | ✗ | contour.pos_.push_back(vec2(A.x, B.y)); | |
| 330 | ✗ | contour.vid_.push_back(int(pts.size())); | |
| 331 | ✗ | pts.push_back(contour.pos_.back()); | |
| 332 | |||
| 333 | ✗ | for (int j = int(B.y - 1); j > int(A.y); j--) { | |
| 334 | ✗ | contour.vid_.push_back(contour.vid(2*singularity_index-contour.pos_.size())); | |
| 335 | ✗ | contour.pos_.push_back(vec2(A.x, double(j))); | |
| 336 | } | ||
| 337 | ✗ | contour.compute_angu(); | |
| 338 | ✗ | contour.angu(singularity_index) = 2 - sing_valence; | |
| 339 | ✗ | return int(singularity_index); | |
| 340 | } | ||
| 341 | |||
| 342 | |||
| 343 | void export_debug_mesh() { | ||
| 344 | Mesh outm; | ||
| 345 | outm.vertices.create_vertices(contour.pos_.size()); | ||
| 346 | vector<index_t> vid(contour.pos_.size()); | ||
| 347 | FOR(i, contour.pos_.size()) { | ||
| 348 | vid[i] = i; | ||
| 349 | X(&outm)[i] = vec3(contour.pos(i)[0], contour.pos(i)[1], 0); | ||
| 350 | } | ||
| 351 | Attribute<int> angu_attr(outm.vertices.attributes(), "angu"); | ||
| 352 | FOR(i, contour.pos_.size()) angu_attr[i] = contour.angu(i); | ||
| 353 | |||
| 354 | outm.facets.create_polygon(vid); | ||
| 355 | mesh_save(outm, "C:/DATA/2dcontours/contour2D"+ String::to_string(export_debug_mesh_id ++) +".geogram"); | ||
| 356 | } | ||
| 357 | |||
| 358 | |||
| 359 | ✗ | bool try_to_punch() { | |
| 360 | ✗ | if (contour.pos_.size() < 4) return false; | |
| 361 | ✗ | FOR(v, contour.pos_.size()) { | |
| 362 | ✗ | if (contour.angu(v) != 1) continue; | |
| 363 | // cut ear | ||
| 364 | ✗ | if (contour.angu(v + 1) == 1) { | |
| 365 | ✗ | FOR(s, 4) quads.push_back(index_t(contour.vid(int(v) - 1 + int(s)))); | |
| 366 | ✗ | contour.remove(int(v)); | |
| 367 | ✗ | contour.remove(int(v)); | |
| 368 | ✗ | contour.angu(int(v) - 1)++; | |
| 369 | ✗ | contour.angu(int(v)) ++; | |
| 370 | ✗ | return true; | |
| 371 | } | ||
| 372 | // add new point | ||
| 373 | ✗ | vec2 npos = contour.pos(v - 1) + contour.pos(v + 1) - contour.pos(v); | |
| 374 | ✗ | bool conflict = false; | |
| 375 | ✗ | FOR(vv, contour.pos_.size()) if (vv != v && (contour.pos(vv) - npos).length2() < .1) conflict = true; | |
| 376 | ✗ | if (conflict) continue; | |
| 377 | ✗ | FOR(s,3) quads.push_back(index_t(contour.vid(int(v) - 1+int(s)))); | |
| 378 | ✗ | quads.push_back(index_t(pts.size())); | |
| 379 | ✗ | contour.vid(v) = int(pts.size()); | |
| 380 | ✗ | pts.push_back(npos); | |
| 381 | ✗ | contour.pos(v) = npos; | |
| 382 | ✗ | contour.angu(v-1) ++; | |
| 383 | ✗ | contour.angu(v) = -1; | |
| 384 | ✗ | contour.angu(v+1) ++; | |
| 385 | ✗ | return true; | |
| 386 | } | ||
| 387 | ✗ | return false; | |
| 388 | } | ||
| 389 | |||
| 390 | |||
| 391 | ✗ | bool apply(vector<int>& angu,int sing_valence) { | |
| 392 | ✗ | int singularity_index=-1; | |
| 393 | ✗ | index_t border_size = pts.size(); | |
| 394 | ✗ | if (sing_valence == 3|| sing_valence == 5) { | |
| 395 | ✗ | singularity_index = init_contour(angu, sing_valence); | |
| 396 | ✗ | } else if (sing_valence == 4) { | |
| 397 | ✗ | contour.resize(pts.size()); | |
| 398 | ✗ | vec2 dir(1, 0); | |
| 399 | ✗ | contour.pos(0) = vec2(0, 0); | |
| 400 | ✗ | for (index_t v = 1; v < pts.size();v++) { | |
| 401 | ✗ | contour.pos(v) = contour.pos(v - 1) + dir; | |
| 402 | ✗ | if (angu[v]< 0) dir = -(R90*dir); | |
| 403 | ✗ | if (angu[v]> 0) dir = R90*dir; | |
| 404 | } | ||
| 405 | ✗ | FOR(v, pts.size()) contour.vid(v) = int(v); | |
| 406 | ✗ | FOR(v, pts.size()) contour.angu(v) = angu[v]; | |
| 407 | ✗ | if ((contour.pos_.back() + dir - contour.pos_.front()).length2() > .1) return false; // check that it is closed | |
| 408 | } | ||
| 409 | else { | ||
| 410 | ✗ | return false; | |
| 411 | } | ||
| 412 | |||
| 413 | //plop("valok"); | ||
| 414 | //plop(sing_valence); | ||
| 415 | //plop(singularity_index); | ||
| 416 | //plop(border_size); | ||
| 417 | //plop(contour.pos_.size()); | ||
| 418 | ✗ | vector<vec2> theta_r(border_size); | |
| 419 | ✗ | vec2 O(.5, .5); | |
| 420 | ✗ | if (singularity_index!=-1) O= contour.pos(singularity_index); | |
| 421 | ✗ | FOR(v, border_size) theta_r[v][1] = (contour.pos(v) - O).length(); | |
| 422 | ✗ | theta_r[0][0] = 0; | |
| 423 | ✗ | FOR(v, border_size - 1) { | |
| 424 | ✗ | theta_r[v + 1][0] = theta_r[v][0] + atan2(det(contour.pos(v) - O, contour.pos(v + 1) - O), dot(contour.pos(v) - O, contour.pos(v + 1) - O)); | |
| 425 | } | ||
| 426 | |||
| 427 | ✗ | FOR(v, border_size) FOR(vv, border_size) if (vv != v && (theta_r[vv] - theta_r[v]).length2() < .0001) { | |
| 428 | ✗ | FOR(w, theta_r.size()) plop(theta_r[w]); | |
| 429 | ✗ | return false; | |
| 430 | } | ||
| 431 | |||
| 432 | |||
| 433 | ✗ | while (try_to_punch()) { | |
| 434 | //export_debug_mesh(); | ||
| 435 | //plop(export_debug_mesh_id); | ||
| 436 | ✗ | if (quads.size() > 1000) geo_assert_not_reached; | |
| 437 | } | ||
| 438 | |||
| 439 | ✗ | nlNewContext(); | |
| 440 | ✗ | nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); | |
| 441 | ✗ | nlSolverParameteri(NL_NB_VARIABLES, NLint(2*pts.size())); | |
| 442 | ✗ | nlBegin(NL_SYSTEM); | |
| 443 | ✗ | FOR(v, border_size) FOR(d,2){ | |
| 444 | ✗ | nlSetVariable(v * 2 +d, pts[v][d]); | |
| 445 | ✗ | nlLockVariable(v * 2 + d); | |
| 446 | } | ||
| 447 | ✗ | nlBegin(NL_MATRIX); | |
| 448 | ✗ | FOR(q, quads.size() / 4) FOR(e, 4) FOR(d, 2) { | |
| 449 | ✗ | nlBegin(NL_ROW); | |
| 450 | ✗ | nlCoefficient(quads[4 * q + e]*2+d, -1.); | |
| 451 | ✗ | nlCoefficient(quads[4 * q + ((e+1)%4)] * 2 + d, 1.); | |
| 452 | ✗ | nlEnd(NL_ROW); | |
| 453 | } | ||
| 454 | |||
| 455 | ✗ | nlEnd(NL_MATRIX); | |
| 456 | ✗ | nlEnd(NL_SYSTEM); | |
| 457 | ✗ | nlSolve(); | |
| 458 | ✗ | FOR(v, pts.size()) FOR(d, 2) pts[v][d]= nlGetVariable(2*v+d); | |
| 459 | ✗ | nlDeleteContext(nlGetCurrent()); | |
| 460 | //export_debug_mesh(); | ||
| 461 | ✗ | if (contour.pos_.size() > 2) { pts.resize(index_t(border_size)); quads.clear(); return false; } | |
| 462 | |||
| 463 | ✗ | if (!Poly2d(pts).quads_are_valid(quads)) { | |
| 464 | //export_debug_mesh(); | ||
| 465 | ✗ | FOR(i, contour.pos_.size()) plop(contour.pos(i)); | |
| 466 | } | ||
| 467 | ✗ | return Poly2d(pts).quads_are_valid(quads); | |
| 468 | |||
| 469 | ✗ | } | |
| 470 | mat2 R90; | ||
| 471 | |||
| 472 | vector<vec2>& pts; | ||
| 473 | vector<index_t>& quads; | ||
| 474 | Contour2D contour; | ||
| 475 | }; | ||
| 476 | |||
| 477 | |||
| 478 | |||
| 479 | |||
| 480 | |||
| 481 | |||
| 482 | |||
| 483 | |||
| 484 | |||
| 485 | |||
| 486 | |||
| 487 | |||
| 488 | |||
| 489 | ✗ | bool Poly2d::try_quad_cover(vector<index_t>& quads) { | |
| 490 | ✗ | vector<int> angu(pts.size(), 1); | |
| 491 | ✗ | int sing_valence = 0; | |
| 492 | ✗ | FOR(v, pts.size()) { | |
| 493 | ✗ | vec2 d0 = aupp(v, pts) - aupp(v - 1, pts); | |
| 494 | ✗ | vec2 d1 = aupp(v + 1, pts) - aupp(v, pts); | |
| 495 | ✗ | double angle = atan2(det(d0, d1), dot(d0, d1)); | |
| 496 | ✗ | if (angle < M_PI / 4.) angu[v] = 0; | |
| 497 | ✗ | if (angle < -M_PI / 4.) angu[v] = -1; | |
| 498 | ✗ | sing_valence += angu[v]; | |
| 499 | } | ||
| 500 | //plop("try_quad_cover"); | ||
| 501 | //dump_contour(); | ||
| 502 | |||
| 503 | ✗ | QuadrangulateWithOneSingularity doit(pts,quads); | |
| 504 | ✗ | if (doit.apply(angu, sing_valence)) return true; | |
| 505 | |||
| 506 | ✗ | return false; | |
| 507 | ✗ | } | |
| 508 | |||
| 509 | |||
| 510 | |||
| 511 | |||
| 512 | |||
| 513 | ✗ | bool Poly2d::try_quadrangulate(vector<index_t>& quads) { | |
| 514 | ✗ | bool verbose = false; | |
| 515 | |||
| 516 | ✗ | index_t nbv = pts.size(); | |
| 517 | ✗ | if (verbose) plop(nbv); | |
| 518 | ✗ | if (nbv < 4) return false; | |
| 519 | ✗ | if (nbv == 4) { | |
| 520 | ✗ | FOR(v, 4) quads.push_back(v); | |
| 521 | ✗ | if (!quads_are_valid(quads)) { GEO::Logger::out("HexDom") << "FAIL" << std::endl; return false; } | |
| 522 | ✗ | return true; | |
| 523 | } | ||
| 524 | |||
| 525 | |||
| 526 | ✗ | if (nbv % 2 != 0) { | |
| 527 | ✗ | GEO::Logger::out("HexDom") << "There is no way to quadrangulate a surface with an odd number of boundary edges" << std::endl; | |
| 528 | ✗ | return false; | |
| 529 | } | ||
| 530 | |||
| 531 | |||
| 532 | ✗ | return try_quad_cover(quads); | |
| 533 | |||
| 534 | /* | ||
| 535 | // precompute a few things | ||
| 536 | vector<double> angle(nbv); | ||
| 537 | vector<double> length(nbv); | ||
| 538 | double ave_length = 0; | ||
| 539 | FOR(i, nbv) { | ||
| 540 | vec2 P[3]; | ||
| 541 | FOR(p, 3) P[p] = aupp(i + p - 1, pts); | ||
| 542 | angle[i] = (180. / M_PI)*atan2(det(P[1] - P[0], P[2] - P[1]), dot(P[1] - P[0], P[2] - P[1])); | ||
| 543 | if (verbose)GEO::Logger::out("HexDom") << "i= " << i << "angle = " << angle[i] << std::endl; | ||
| 544 | length[i] = (P[1] - P[0]).length() / double(nbv); | ||
| 545 | ave_length += length[i]; | ||
| 546 | } | ||
| 547 | plop("gna"); | ||
| 548 | |||
| 549 | |||
| 550 | // define outputs of the search | ||
| 551 | index_t start = index_t(-1); | ||
| 552 | index_t end = index_t(-1); | ||
| 553 | index_t nb_nv_pts = index_t(-1); | ||
| 554 | double best_score = 0; | ||
| 555 | |||
| 556 | |||
| 557 | |||
| 558 | index_t dec = parity_of_original_points(); | ||
| 559 | plop("gna"); | ||
| 560 | |||
| 561 | FOR(test_start, nbv) { | ||
| 562 | plop(test_start); | ||
| 563 | index_t test_end; | ||
| 564 | index_t test_nb_nv_pts; | ||
| 565 | double test_score; | ||
| 566 | |||
| 567 | plop("gna"); | ||
| 568 | |||
| 569 | FOR(d, nbv - 5) { | ||
| 570 | plop(d); | ||
| 571 | |||
| 572 | test_end = test_start + d + 3; | ||
| 573 | |||
| 574 | vec2 A[3]; | ||
| 575 | FOR(i, 3) A[i] = aupp(int(test_start + i) - 1, pts); | ||
| 576 | vec2 B[3]; | ||
| 577 | FOR(i, 3) B[i] = aupp(int(test_end + i) - 1, pts); | ||
| 578 | vec2 nA1A2 = normalize(A[2] - A[1]); | ||
| 579 | vec2 nA1A0 = normalize(A[0] - A[1]); | ||
| 580 | vec2 nB1B2 = normalize(B[2] - B[1]); | ||
| 581 | vec2 nB1B0 = normalize(B[0] - B[1]); | ||
| 582 | vec2 nAB = normalize(B[1] - A[1]); | ||
| 583 | vec2 nBA = -nAB; | ||
| 584 | |||
| 585 | double worst_det = 1; | ||
| 586 | worst_det = std::min(worst_det, det(nA1A2, nAB)); | ||
| 587 | worst_det = std::min(worst_det, det(nAB, nA1A0)); | ||
| 588 | worst_det = std::min(worst_det, det(nB1B2, nBA)); | ||
| 589 | worst_det = std::min(worst_det, det(nBA, nB1B0)); | ||
| 590 | |||
| 591 | test_score = worst_det; | ||
| 592 | |||
| 593 | |||
| 594 | double AB_relative_length = floor((B[1] - A[1]).length() / ave_length); | ||
| 595 | test_nb_nv_pts = index_t(std::max(0, int(AB_relative_length) - 1)); | ||
| 596 | |||
| 597 | if (test_nb_nv_pts % 2 != int(d % 2)) { | ||
| 598 | if (test_nb_nv_pts == 0) test_nb_nv_pts = 1; else test_nb_nv_pts--; | ||
| 599 | } | ||
| 600 | |||
| 601 | if (angle[test_start] < 1) test_score += 1; | ||
| 602 | if (angle[test_end] < 1) test_score += 1; | ||
| 603 | if (angle[test_start] < -45) test_score += 2; | ||
| 604 | if (angle[test_end] < -45) test_score += 2; | ||
| 605 | if ((test_start % 2) == 1 - dec) test_score -= 10; | ||
| 606 | if ((test_end % 2) == 1 - dec) test_score -= 10; | ||
| 607 | test_nb_nv_pts = 1; | ||
| 608 | |||
| 609 | if (best_score < test_score) { | ||
| 610 | bool can_cut = true; | ||
| 611 | FOR(dd, nbv) { | ||
| 612 | index_t ind = test_start + dd; | ||
| 613 | if (ind > test_start && ind < test_end) | ||
| 614 | can_cut = can_cut && det(nAB, aupp(ind, pts) - A[1]) < 0; | ||
| 615 | if (ind > test_end) | ||
| 616 | can_cut = can_cut && det(nAB, aupp(ind, pts) - A[1]) > 0; | ||
| 617 | } | ||
| 618 | if (verbose) | ||
| 619 | std::cerr << "can_cut = " << can_cut << " test_score = " << test_score | ||
| 620 | << " test_start = " << test_start << " test_end = " << test_end | ||
| 621 | << " test_nb_nv_pts = " << test_nb_nv_pts << std::endl; | ||
| 622 | if (can_cut) { | ||
| 623 | start = test_start; | ||
| 624 | end = test_end; | ||
| 625 | nb_nv_pts = test_nb_nv_pts; | ||
| 626 | best_score = test_score; | ||
| 627 | } | ||
| 628 | } | ||
| 629 | } | ||
| 630 | |||
| 631 | |||
| 632 | } | ||
| 633 | |||
| 634 | plop("gna"); | ||
| 635 | |||
| 636 | if (nbv > 8) | ||
| 637 | if (nb_nv_pts != index_t(-1)) { | ||
| 638 | if (verbose)GEO::Logger::out("HexDom") << "remove quad strip from " << start << " with score = " << best_score << " with nbpts" << nb_nv_pts << std::endl; | ||
| 639 | |||
| 640 | vector<index_t> global_vid[2]; // gives indices in "pts" from indices in "poly[i]" | ||
| 641 | |||
| 642 | // fill both half with existing points | ||
| 643 | //int end = start + nb_nv_pts + 3; | ||
| 644 | FOR(d, end - start + 1) global_vid[0].push_back((start + d) % nbv); | ||
| 645 | FOR(d, nbv - (end - start) + 1) global_vid[1].push_back((end + d) % nbv); | ||
| 646 | |||
| 647 | |||
| 648 | // add new vertices along the cut | ||
| 649 | FOR(i, nb_nv_pts) global_vid[0].push_back(nbv + i); | ||
| 650 | FOR(i, nb_nv_pts) global_vid[1].push_back(nbv + (nb_nv_pts - 1 - i)); | ||
| 651 | FOR(i, nb_nv_pts) { | ||
| 652 | double c = 1.0 - double(i + 1) / double(nb_nv_pts + 1); | ||
| 653 | pts.push_back((1. - c)*pts[start] + c*pts[end% nbv]); | ||
| 654 | } | ||
| 655 | |||
| 656 | // solve on two halves | ||
| 657 | vector<vec2> poly[2]; | ||
| 658 | FOR(i, 2) FOR(fv, global_vid[i].size()) poly[i].push_back(pts[global_vid[i][fv]]); | ||
| 659 | |||
| 660 | vector<index_t> poly_quad[2]; | ||
| 661 | FOR(i, 2) if (!Poly2d(poly[i]).try_quadrangulate(poly_quad[i])) return false; | ||
| 662 | // add new pts to global | ||
| 663 | FOR(i, 2) for (index_t d = global_vid[i].size(); d < poly[i].size(); d++) { | ||
| 664 | global_vid[i].push_back(pts.size()); | ||
| 665 | pts.push_back(poly[i][d]); | ||
| 666 | } | ||
| 667 | FOR(i, 2) FOR(qu, poly_quad[i].size()) quads.push_back(global_vid[i][poly_quad[i][qu]]); | ||
| 668 | |||
| 669 | if (!quads_are_valid(quads)) { GEO::Logger::out("HexDom") << "FAIL remove quad strip" << std::endl; return false; } | ||
| 670 | return true; | ||
| 671 | } | ||
| 672 | |||
| 673 | |||
| 674 | plop("gna"); | ||
| 675 | |||
| 676 | |||
| 677 | if (verbose) GEO::Logger::out("HexDom") << "middle_point_quadrangulate(quads)" << std::endl; | ||
| 678 | middle_point_quadrangulate(quads); | ||
| 679 | if (!quads_are_valid(quads)) { GEO::Logger::out("HexDom") << "FAIL middle_point_quadrangulate" << std::endl; return false; } | ||
| 680 | |||
| 681 | return true; | ||
| 682 | */ | ||
| 683 | } | ||
| 684 | |||
| 685 | /*****************************************************************************************************/ | ||
| 686 | |||
| 687 | ✗ | vec3 Poly3d::barycenter() { | |
| 688 | ✗ | vec3 bary(0, 0, 0); | |
| 689 | ✗ | FOR(fv, pts.size()) { | |
| 690 | ✗ | bary = bary + (1. / double(pts.size()))*pts[fv]; | |
| 691 | } | ||
| 692 | ✗ | return bary; | |
| 693 | } | ||
| 694 | |||
| 695 | ✗ | vec3 Poly3d::normal() { | |
| 696 | ✗ | vec3 n(0, 0, 0); | |
| 697 | ✗ | vec3 bary = barycenter(); | |
| 698 | ✗ | FOR(fv, pts.size()) { | |
| 699 | ✗ | n = n + cross(pts[fv] - bary, pts[next_mod(fv, pts.size())] - bary); | |
| 700 | // plop(n); | ||
| 701 | } | ||
| 702 | ✗ | n = normalize(n); | |
| 703 | ✗ | return n; | |
| 704 | } | ||
| 705 | |||
| 706 | |||
| 707 | ✗ | bool Poly3d::try_triangulate_minweight(vector<index_t>& triangles) { | |
| 708 | ✗ | index_t nbv = pts.size(); | |
| 709 | ✗ | if (nbv == 3) { | |
| 710 | ✗ | FOR(v, 3) { | |
| 711 | ✗ | triangles.push_back(v); | |
| 712 | } | ||
| 713 | ✗ | return true; | |
| 714 | } | ||
| 715 | ✗ | geo_assert(nbv > 3); | |
| 716 | |||
| 717 | ✗ | vector<vec2> pts2d; | |
| 718 | ✗ | Basis3d b(normal()); | |
| 719 | ✗ | FOR(fv, nbv) { | |
| 720 | ✗ | pts2d.push_back(b.project_xy(pts[fv])); | |
| 721 | } | ||
| 722 | |||
| 723 | ✗ | return Poly2d(pts2d).try_triangulate_minweight(triangles); | |
| 724 | ✗ | } | |
| 725 | |||
| 726 | /** | ||
| 727 | * WARNING: it may introduce new vertices in pts | ||
| 728 | */ | ||
| 729 | ✗ | bool Poly3d::try_quadrangulate(vector<index_t>& quads) { | |
| 730 | ✗ | index_t nbv = pts.size(); | |
| 731 | ✗ | if (nbv < 4) return false; | |
| 732 | ✗ | vec3 G = barycenter(); | |
| 733 | ✗ | if (normal().length2() < 1e-20) return false; | |
| 734 | ✗ | Basis3d b(normal()); | |
| 735 | |||
| 736 | ✗ | vector<vec2> pts2d; | |
| 737 | ✗ | FOR(fv, nbv) pts2d.push_back(b.project_xy(pts[fv] - G)); | |
| 738 | ✗ | Poly2d p2d(pts2d); | |
| 739 | ✗ | if (!p2d.try_quadrangulate(quads)) { | |
| 740 | //dump_contour(); | ||
| 741 | ✗ | return false; | |
| 742 | } | ||
| 743 | ✗ | for (index_t i = pts.size(); i < p2d.pts.size(); i++) | |
| 744 | ✗ | pts.push_back(G + b.un_project_xy(p2d.pts[i])); | |
| 745 | ✗ | return true; | |
| 746 | ✗ | } | |
| 747 | |||
| 748 | } | ||
| 749 |