| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2022 Inria | ||
| 3 | * All rights reserved. | ||
| 4 | * | ||
| 5 | * Redistribution and use in source and binary forms, with or without | ||
| 6 | * modification, are permitted provided that the following conditions are met: | ||
| 7 | * | ||
| 8 | * * Redistributions of source code must retain the above copyright notice, | ||
| 9 | * this list of conditions and the following disclaimer. | ||
| 10 | * * Redistributions in binary form must reproduce the above copyright notice, | ||
| 11 | * this list of conditions and the following disclaimer in the documentation | ||
| 12 | * and/or other materials provided with the distribution. | ||
| 13 | * * Neither the name of the ALICE Project-Team nor the names of its | ||
| 14 | * contributors may be used to endorse or promote products derived from this | ||
| 15 | * software without specific prior written permission. | ||
| 16 | * | ||
| 17 | * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | ||
| 18 | * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
| 19 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
| 20 | * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE | ||
| 21 | * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | ||
| 22 | * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | ||
| 23 | * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | ||
| 24 | * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN | ||
| 25 | * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | ||
| 26 | * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE | ||
| 27 | * POSSIBILITY OF SUCH DAMAGE. | ||
| 28 | * | ||
| 29 | * Contact: Bruno Levy | ||
| 30 | * | ||
| 31 | * https://www.inria.fr/fr/bruno-levy | ||
| 32 | * | ||
| 33 | * Inria, | ||
| 34 | * Domaine de Voluceau, | ||
| 35 | * 78150 Le Chesnay - Rocquencourt | ||
| 36 | * FRANCE | ||
| 37 | * | ||
| 38 | */ | ||
| 39 | |||
| 40 | #include <geogram/parameterization/mesh_LSCM.h> | ||
| 41 | #include <geogram/mesh/mesh.h> | ||
| 42 | #include <geogram/mesh/mesh_geometry.h> | ||
| 43 | #include <geogram/bibliography/bibliography.h> | ||
| 44 | #include <geogram/NL/nl.h> | ||
| 45 | #include <geogram/NL/nl_ext.h> | ||
| 46 | #include <algorithm> | ||
| 47 | |||
| 48 | namespace { | ||
| 49 | using namespace GEO; | ||
| 50 | |||
| 51 | |||
| 52 | /** | ||
| 53 | * \brief Computes Least Squares Conformal Maps in least squares or | ||
| 54 | * spectral mode. | ||
| 55 | * \details The method is described in the following references: | ||
| 56 | * - Least Squares Conformal Maps, Levy, Petitjean, Ray, Maillot, ACM | ||
| 57 | * SIGGRAPH, 2002 | ||
| 58 | * - Spectral Conformal Parameterization, Mullen, Tong, Alliez, Desbrun, | ||
| 59 | * Computer Graphics Forum (SGP conf. proc.), 2008 | ||
| 60 | */ | ||
| 61 | class LSCM { | ||
| 62 | public: | ||
| 63 | |||
| 64 | /** | ||
| 65 | * \brief LSCM constructor | ||
| 66 | * \param[in] M a reference to the mesh. It needs to correspond to a | ||
| 67 | * topological disk (open surface with one border and no handle). | ||
| 68 | * \param[in] tex_coord a vector property of dimension 2 where to | ||
| 69 | * store the texture coordinates. | ||
| 70 | */ | ||
| 71 | ✗ | LSCM(Mesh& M, Attribute<double>& tex_coord, Attribute<double>& angle) : | |
| 72 | ✗ | mesh_(M), tex_coord_(tex_coord), angle_(angle), eigen_(0) { | |
| 73 | ✗ | geo_assert(tex_coord.dimension() == 2); | |
| 74 | ✗ | locked_1_ = NO_INDEX; | |
| 75 | ✗ | locked_2_ = NO_INDEX; | |
| 76 | ✗ | verbose_ = false; | |
| 77 | ✗ | } | |
| 78 | |||
| 79 | /** | ||
| 80 | * \brief Enables or disables messages. | ||
| 81 | * \param[in] x if true, messages are displayed on the console | ||
| 82 | * with statistics. Default is non-verbose. | ||
| 83 | */ | ||
| 84 | ✗ | void set_verbose(bool x) { | |
| 85 | ✗ | verbose_ = x; | |
| 86 | ✗ | } | |
| 87 | |||
| 88 | /** | ||
| 89 | * \brief Sets whether spectral mode is used. | ||
| 90 | * \details In default mode, the trivial solution (all vertices to zero) | ||
| 91 | * is avoided by locking two vertices (that are as "extremal" | ||
| 92 | * as possible). In spectral mode, the trivial solution is avoided by | ||
| 93 | * finding the first minimizer that is orthogonal to it (more elegant, | ||
| 94 | * but more costly). | ||
| 95 | */ | ||
| 96 | ✗ | void set_spectral(bool x) { | |
| 97 | ✗ | spectral_ = x; | |
| 98 | ✗ | } | |
| 99 | |||
| 100 | /** | ||
| 101 | * \brief Computes the least squares conformal map and stores it in | ||
| 102 | * the texture coordinates of the mesh. | ||
| 103 | * \details Outline of the algorithm (steps 1,2,3 are not used | ||
| 104 | * in spectral mode): | ||
| 105 | * - 1) Find an initial solution by projecting on a plane | ||
| 106 | * - 2) Lock two vertices of the mesh | ||
| 107 | * - 3) Copy the initial u,v coordinates to OpenNL | ||
| 108 | * - 4) Construct the LSCM equation with OpenNL | ||
| 109 | * - 5) Solve the equation with OpenNL | ||
| 110 | * - 6) Copy OpenNL solution to the u,v coordinates | ||
| 111 | */ | ||
| 112 | |||
| 113 | ✗ | void apply() { | |
| 114 | |||
| 115 | ✗ | geo_cite("DBLP:journals/tog/LevyPRM02"); | |
| 116 | ✗ | if(spectral_) { | |
| 117 | ✗ | geo_cite("DBLP:journals/cgf/MullenTAD08"); | |
| 118 | } | ||
| 119 | |||
| 120 | ✗ | const int nb_eigens = 10; | |
| 121 | ✗ | nlNewContext(); | |
| 122 | ✗ | NLuint nb_vertices = NLuint(mesh_.vertices.nb()); | |
| 123 | |||
| 124 | ✗ | if(spectral_) { | |
| 125 | ✗ | if(nlInitExtension("ARPACK")) { | |
| 126 | ✗ | if(verbose_) { | |
| 127 | ✗ | Logger::out("LSCM") << "ARPACK extension initialized" | |
| 128 | ✗ | << std::endl; | |
| 129 | } | ||
| 130 | ✗ | nlEigenSolverParameteri(NL_EIGEN_SOLVER, NL_ARPACK_EXT); | |
| 131 | ✗ | nlEigenSolverParameteri(NL_NB_EIGENS, nb_eigens); | |
| 132 | ✗ | if(verbose_) { | |
| 133 | ✗ | nlEnable(NL_VERBOSE); | |
| 134 | } | ||
| 135 | } else { | ||
| 136 | ✗ | if(verbose_) { | |
| 137 | ✗ | Logger::out("LSCM") | |
| 138 | ✗ | << "Could not initialize ARPACK extension" | |
| 139 | ✗ | << std::endl; | |
| 140 | ✗ | Logger::out("LSCM") | |
| 141 | ✗ | << "Falling back to least squares mode" | |
| 142 | ✗ | << std::endl; | |
| 143 | } | ||
| 144 | ✗ | spectral_ = false; | |
| 145 | } | ||
| 146 | } else { | ||
| 147 | #ifndef GEO_OS_ANDROID | ||
| 148 | ✗ | if( | |
| 149 | ✗ | nb_vertices <= 200000 && ( | |
| 150 | ✗ | nlExtensionIsInitialized("SUPERLU") || | |
| 151 | ✗ | nlInitExtension("SUPERLU") | |
| 152 | ) | ||
| 153 | ) { | ||
| 154 | ✗ | if(verbose_) { | |
| 155 | ✗ | Logger::out("LSCM") << "using SUPERLU" | |
| 156 | ✗ | << std::endl; | |
| 157 | } | ||
| 158 | ✗ | nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT); | |
| 159 | ✗ | } else if( | |
| 160 | ✗ | nlExtensionIsInitialized("AMGCL") || | |
| 161 | ✗ | nlInitExtension("AMGCL") | |
| 162 | ) { | ||
| 163 | ✗ | if(verbose_) { | |
| 164 | ✗ | nlEnable(NL_VERBOSE); | |
| 165 | ✗ | Logger::out("LSCM") << "using AMGCL" | |
| 166 | ✗ | << std::endl; | |
| 167 | } | ||
| 168 | ✗ | nlSolverParameteri(NL_SOLVER, NL_AMGCL_EXT); | |
| 169 | } else | ||
| 170 | #endif | ||
| 171 | { | ||
| 172 | ✗ | if(verbose_) { | |
| 173 | ✗ | nlEnable(NL_VERBOSE); | |
| 174 | ✗ | Logger::out("LSCM") << "using JacobiCG" | |
| 175 | ✗ | << std::endl; | |
| 176 | } | ||
| 177 | } | ||
| 178 | } | ||
| 179 | ✗ | if(!spectral_) { | |
| 180 | ✗ | project(); | |
| 181 | } | ||
| 182 | ✗ | nlSolverParameteri(NL_NB_VARIABLES, NLint(2*nb_vertices)); | |
| 183 | ✗ | nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); | |
| 184 | ✗ | nlSolverParameteri(NL_MAX_ITERATIONS, NLint(5*nb_vertices)); | |
| 185 | ✗ | if(spectral_) { | |
| 186 | ✗ | nlSolverParameterd(NL_THRESHOLD, 0.0); | |
| 187 | } else { | ||
| 188 | ✗ | nlSolverParameterd(NL_THRESHOLD, 1e-10); | |
| 189 | } | ||
| 190 | ✗ | nlBegin(NL_SYSTEM); | |
| 191 | ✗ | mesh_to_solver(); | |
| 192 | ✗ | nlBegin(NL_MATRIX); | |
| 193 | ✗ | setup_lscm(); | |
| 194 | ✗ | nlEnd(NL_MATRIX); | |
| 195 | ✗ | nlEnd(NL_SYSTEM); | |
| 196 | ✗ | if(verbose_) { | |
| 197 | ✗ | Logger::out("LSCM") << "Solving ..." << std::endl; | |
| 198 | } | ||
| 199 | |||
| 200 | ✗ | if(spectral_) { | |
| 201 | ✗ | nlEigenSolve(); | |
| 202 | ✗ | if(verbose_) { | |
| 203 | ✗ | for(NLuint i=0; i<nb_eigens; ++i) { | |
| 204 | ✗ | Logger::out("LSCM") << "[" << i << "] " | |
| 205 | ✗ | << nlGetEigenValue(i) << std::endl; | |
| 206 | } | ||
| 207 | } | ||
| 208 | |||
| 209 | // Find first "non-zero" eigenvalue | ||
| 210 | ✗ | double small_eigen = ::fabs(nlGetEigenValue(0)) ; | |
| 211 | ✗ | eigen_ = 1; | |
| 212 | ✗ | for(NLuint i=1; i<nb_eigens; ++i) { | |
| 213 | ✗ | if(::fabs(nlGetEigenValue(i)) / small_eigen > 1e3) { | |
| 214 | ✗ | eigen_ = i ; | |
| 215 | ✗ | break ; | |
| 216 | } | ||
| 217 | } | ||
| 218 | } else{ | ||
| 219 | ✗ | nlSolve(); | |
| 220 | } | ||
| 221 | |||
| 222 | ✗ | solver_to_mesh(); | |
| 223 | ✗ | normalize_uv(); | |
| 224 | |||
| 225 | ✗ | if(!spectral_) { | |
| 226 | ✗ | if(verbose_) { | |
| 227 | double time; | ||
| 228 | NLint iterations; | ||
| 229 | ✗ | nlGetDoublev(NL_ELAPSED_TIME, &time); | |
| 230 | ✗ | nlGetIntegerv(NL_USED_ITERATIONS, &iterations); | |
| 231 | ✗ | Logger::out("LSCM") << "Solver time: " << time << std::endl; | |
| 232 | ✗ | Logger::out("LSCM") << "Used iterations: " | |
| 233 | ✗ | << iterations << std::endl; | |
| 234 | } | ||
| 235 | } | ||
| 236 | ✗ | nlDeleteContext(nlGetCurrent()); | |
| 237 | ✗ | } | |
| 238 | |||
| 239 | protected: | ||
| 240 | |||
| 241 | /** | ||
| 242 | * \brief Creates the LSCM equations in OpenNL. | ||
| 243 | */ | ||
| 244 | ✗ | void setup_lscm() { | |
| 245 | ✗ | for(NLuint f=0; f<mesh_.facets.nb(); ++f) { | |
| 246 | ✗ | setup_lscm(f); | |
| 247 | } | ||
| 248 | ✗ | } | |
| 249 | |||
| 250 | /** | ||
| 251 | * \brief Creates the LSCM equations in OpenNL, related | ||
| 252 | * with a given facet. | ||
| 253 | * \param[in] f the index of the facet. | ||
| 254 | * \details no-need to triangulate the facet, | ||
| 255 | * we do that "virtually", by creating triangles | ||
| 256 | * radiating around vertex 0 of the facet. | ||
| 257 | * (however, this may be invalid for concave facets) | ||
| 258 | */ | ||
| 259 | ✗ | void setup_lscm(NLuint f) { | |
| 260 | ✗ | NLuint nv = NLuint(mesh_.facets.nb_vertices(f)); | |
| 261 | ✗ | if(angle_.is_bound()) { | |
| 262 | ✗ | index_t c0 = mesh_.facets.corners_begin(f); | |
| 263 | ✗ | double a0 = angle_[c0]; | |
| 264 | ✗ | for(NLuint i=1; i<nv-1; ++i) { | |
| 265 | ✗ | double ai = angle_[c0+i]; | |
| 266 | ✗ | double aip1 = angle_[c0+1+1]; | |
| 267 | ✗ | setup_conformal_map_relations( | |
| 268 | ✗ | NLuint(mesh_.facets.vertex(f,0)), | |
| 269 | ✗ | NLuint(mesh_.facets.vertex(f,i)), | |
| 270 | ✗ | NLuint(mesh_.facets.vertex(f,i+1)), | |
| 271 | a0, ai, aip1 | ||
| 272 | ); | ||
| 273 | } | ||
| 274 | } else { | ||
| 275 | ✗ | for(NLuint i=1; i<nv-1; ++i) { | |
| 276 | ✗ | setup_conformal_map_relations( | |
| 277 | ✗ | NLuint(mesh_.facets.vertex(f,0)), | |
| 278 | ✗ | NLuint(mesh_.facets.vertex(f,i)), | |
| 279 | ✗ | NLuint(mesh_.facets.vertex(f,i+1)) | |
| 280 | ); | ||
| 281 | } | ||
| 282 | } | ||
| 283 | ✗ | } | |
| 284 | |||
| 285 | /** | ||
| 286 | * \brief Computes the coordinates of the vertices of a triangle | ||
| 287 | * in a local 2D orthonormal basis of the triangle's plane. | ||
| 288 | * \param[in] p0 , p1 , p2 the 3D coordinates of the vertices of | ||
| 289 | * the triangle | ||
| 290 | * \param[out] z0 , z1 , z2 the 2D coordinates of the vertices of | ||
| 291 | * the triangle | ||
| 292 | */ | ||
| 293 | ✗ | static void project_triangle( | |
| 294 | const vec3& p0, | ||
| 295 | const vec3& p1, | ||
| 296 | const vec3& p2, | ||
| 297 | vec2& z0, | ||
| 298 | vec2& z1, | ||
| 299 | vec2& z2 | ||
| 300 | ) { | ||
| 301 | ✗ | vec3 X = p1 - p0; | |
| 302 | ✗ | X = normalize(X); | |
| 303 | ✗ | vec3 Z = cross(X,(p2 - p0)); | |
| 304 | ✗ | Z = normalize(Z); | |
| 305 | ✗ | vec3 Y = cross(Z,X); | |
| 306 | ✗ | const vec3& O = p0; | |
| 307 | |||
| 308 | ✗ | double x0 = 0; | |
| 309 | ✗ | double y0 = 0; | |
| 310 | ✗ | double x1 = (p1 - O).length(); | |
| 311 | ✗ | double y1 = 0; | |
| 312 | ✗ | double x2 = dot((p2 - O),X); | |
| 313 | ✗ | double y2 = dot((p2 - O),Y); | |
| 314 | |||
| 315 | ✗ | z0 = vec2(x0,y0); | |
| 316 | ✗ | z1 = vec2(x1,y1); | |
| 317 | ✗ | z2 = vec2(x2,y2); | |
| 318 | ✗ | } | |
| 319 | |||
| 320 | /** | ||
| 321 | * \brief Creates the LSCM equation in OpenNL, related with | ||
| 322 | * a given triangle, specified by vertex indices. | ||
| 323 | * \param[in] v0 , v1 , v2 the indices of the three vertices of | ||
| 324 | * the triangle. | ||
| 325 | * \details Uses the geometric form of LSCM equation: | ||
| 326 | * (Z1 - Z0)(U2 - U0) = (Z2 - Z0)(U1 - U0) | ||
| 327 | * Where Uk = uk + i.vk is the complex number | ||
| 328 | * corresponding to (u,v) coords | ||
| 329 | * Zk = xk + i.yk is the complex number | ||
| 330 | * corresponding to local (x,y) coords | ||
| 331 | * There is no divide with this expression, | ||
| 332 | * this makes it more numerically stable in | ||
| 333 | * the presence of degenerate triangles. | ||
| 334 | */ | ||
| 335 | ✗ | void setup_conformal_map_relations( | |
| 336 | NLuint v0, NLuint v1, NLuint v2 | ||
| 337 | ) { | ||
| 338 | |||
| 339 | ✗ | const vec3& p0 = mesh_.vertices.point(v0); | |
| 340 | ✗ | const vec3& p1 = mesh_.vertices.point(v1); | |
| 341 | ✗ | const vec3& p2 = mesh_.vertices.point(v2); | |
| 342 | |||
| 343 | ✗ | vec2 z0,z1,z2; | |
| 344 | ✗ | project_triangle(p0,p1,p2,z0,z1,z2); | |
| 345 | ✗ | vec2 z01 = z1 - z0; | |
| 346 | ✗ | vec2 z02 = z2 - z0; | |
| 347 | ✗ | double a = z01.x; | |
| 348 | ✗ | double b = z01.y; | |
| 349 | ✗ | double c = z02.x; | |
| 350 | ✗ | double d = z02.y; | |
| 351 | ✗ | geo_assert(b == 0.0); | |
| 352 | |||
| 353 | // Note : 2*id + 0 --> u | ||
| 354 | // 2*id + 1 --> v | ||
| 355 | ✗ | NLuint u0_id = 2*v0 ; | |
| 356 | ✗ | NLuint v0_id = 2*v0 + 1; | |
| 357 | ✗ | NLuint u1_id = 2*v1 ; | |
| 358 | ✗ | NLuint v1_id = 2*v1 + 1; | |
| 359 | ✗ | NLuint u2_id = 2*v2 ; | |
| 360 | ✗ | NLuint v2_id = 2*v2 + 1; | |
| 361 | |||
| 362 | // Note : rhs = 0 | ||
| 363 | |||
| 364 | // Real part | ||
| 365 | ✗ | nlBegin(NL_ROW); | |
| 366 | ✗ | nlCoefficient(u0_id, -a+c) ; | |
| 367 | ✗ | nlCoefficient(v0_id, b-d) ; | |
| 368 | ✗ | nlCoefficient(u1_id, -c) ; | |
| 369 | ✗ | nlCoefficient(v1_id, d) ; | |
| 370 | ✗ | nlCoefficient(u2_id, a); | |
| 371 | ✗ | nlEnd(NL_ROW); | |
| 372 | |||
| 373 | // Imaginary part | ||
| 374 | ✗ | nlBegin(NL_ROW); | |
| 375 | ✗ | nlCoefficient(u0_id, -b+d); | |
| 376 | ✗ | nlCoefficient(v0_id, -a+c); | |
| 377 | ✗ | nlCoefficient(u1_id, -d); | |
| 378 | ✗ | nlCoefficient(v1_id, -c); | |
| 379 | ✗ | nlCoefficient(v2_id, a); | |
| 380 | ✗ | nlEnd(NL_ROW); | |
| 381 | ✗ | } | |
| 382 | |||
| 383 | /** | ||
| 384 | * \brief Creates the LSCM equation in OpenNL, related with | ||
| 385 | * a given triangle, specified by vertex indices, and with | ||
| 386 | * specified desired angles. | ||
| 387 | * \details This version is used to recoved the u,v coordinates | ||
| 388 | * from the angles computed by ABF++. | ||
| 389 | * \param[in] v0 , v1 , v2 the indices of the three vertices of | ||
| 390 | * the triangle. | ||
| 391 | * \param[in] alpha0 , alpha1 , alpha2 the desired angles at the | ||
| 392 | * three vertices of the triangle | ||
| 393 | */ | ||
| 394 | ✗ | void setup_conformal_map_relations( | |
| 395 | NLuint v0, NLuint v1, NLuint v2, | ||
| 396 | double alpha0, double alpha1, double alpha2 | ||
| 397 | ) { | ||
| 398 | ✗ | const vec3& p0 = mesh_.vertices.point(v0); | |
| 399 | ✗ | const vec3& p1 = mesh_.vertices.point(v1); | |
| 400 | ✗ | const vec3& p2 = mesh_.vertices.point(v2); | |
| 401 | |||
| 402 | ✗ | double scaling = ::sin(alpha1) / ::sin(alpha2) ; | |
| 403 | ✗ | double a = scaling * ::cos(alpha0); | |
| 404 | ✗ | double b = scaling * ::sin(alpha0); | |
| 405 | |||
| 406 | ✗ | double area = Geom::triangle_area(p0,p1,p2) ; | |
| 407 | ✗ | double s = ::sqrt(area) ; | |
| 408 | |||
| 409 | // Note : 2*id + 0 --> u | ||
| 410 | // 2*id + 1 --> v | ||
| 411 | ✗ | NLuint u0_id = 2*v0 ; | |
| 412 | ✗ | NLuint v0_id = 2*v0 + 1; | |
| 413 | ✗ | NLuint u1_id = 2*v1 ; | |
| 414 | ✗ | NLuint v1_id = 2*v1 + 1; | |
| 415 | ✗ | NLuint u2_id = 2*v2 ; | |
| 416 | ✗ | NLuint v2_id = 2*v2 + 1; | |
| 417 | |||
| 418 | // Note : rhs = 0 | ||
| 419 | |||
| 420 | // Real part | ||
| 421 | ✗ | nlRowScaling(s); | |
| 422 | ✗ | nlBegin(NL_ROW); | |
| 423 | ✗ | nlCoefficient(u0_id, 1.0 - a) ; | |
| 424 | ✗ | nlCoefficient(v0_id, b) ; | |
| 425 | ✗ | nlCoefficient(u1_id, a) ; | |
| 426 | ✗ | nlCoefficient(v1_id, -b) ; | |
| 427 | ✗ | nlCoefficient(u2_id, -1.0); | |
| 428 | ✗ | nlEnd(NL_ROW); | |
| 429 | |||
| 430 | // Imaginary part | ||
| 431 | ✗ | nlRowScaling(s); | |
| 432 | ✗ | nlBegin(NL_ROW); | |
| 433 | ✗ | nlCoefficient(u0_id, -b); | |
| 434 | ✗ | nlCoefficient(v0_id, 1.0-a); | |
| 435 | ✗ | nlCoefficient(u1_id, b); | |
| 436 | ✗ | nlCoefficient(v1_id, a); | |
| 437 | ✗ | nlCoefficient(v2_id, -1.0); | |
| 438 | ✗ | nlEnd(NL_ROW); | |
| 439 | ✗ | } | |
| 440 | |||
| 441 | /** | ||
| 442 | * \brief Copies u,v coordinates from OpenNL solver to the mesh. | ||
| 443 | */ | ||
| 444 | ✗ | void solver_to_mesh() { | |
| 445 | ✗ | for(index_t i: mesh_.vertices) { | |
| 446 | ✗ | double u = spectral_ ? nlMultiGetVariable(NLuint(2*i),eigen_) | |
| 447 | ✗ | : nlGetVariable(2*i); | |
| 448 | ✗ | double v = spectral_ ? nlMultiGetVariable(NLuint(2*i+1),eigen_) | |
| 449 | ✗ | : nlGetVariable(2*i+1); | |
| 450 | ✗ | tex_coord_[2*i] = u; | |
| 451 | ✗ | tex_coord_[2*i+1] = v; | |
| 452 | } | ||
| 453 | ✗ | } | |
| 454 | |||
| 455 | /** | ||
| 456 | * \brief Translates and scales tex coords in such a way that they fit | ||
| 457 | * within the unit square. | ||
| 458 | */ | ||
| 459 | ✗ | void normalize_uv() { | |
| 460 | ✗ | double u_min= Numeric::max_float64(); | |
| 461 | ✗ | double v_min= Numeric::max_float64(); | |
| 462 | ✗ | double u_max= -Numeric::max_float64(); | |
| 463 | ✗ | double v_max= -Numeric::max_float64(); | |
| 464 | ✗ | for(NLuint i=0; i<mesh_.vertices.nb(); ++i) { | |
| 465 | ✗ | double u = tex_coord_[2*i]; | |
| 466 | ✗ | double v = tex_coord_[2*i+1]; | |
| 467 | ✗ | u_min = std::min(u_min, u); | |
| 468 | ✗ | v_min = std::min(v_min, v); | |
| 469 | ✗ | u_max = std::max(u_max, u); | |
| 470 | ✗ | v_max = std::max(v_max, v); | |
| 471 | } | ||
| 472 | ✗ | double l = std::max(u_max-u_min,v_max-v_min); | |
| 473 | ✗ | for(NLuint i=0; i<mesh_.vertices.nb(); ++i) { | |
| 474 | ✗ | tex_coord_[2*i] -= u_min; | |
| 475 | ✗ | tex_coord_[2*i] /= l; | |
| 476 | ✗ | tex_coord_[2*i+1] -= v_min; | |
| 477 | ✗ | tex_coord_[2*i+1] /= l; | |
| 478 | } | ||
| 479 | ✗ | } | |
| 480 | |||
| 481 | /** | ||
| 482 | * \brief Tests whether a vertex is locked. | ||
| 483 | * \param[in] v the index of the vertex | ||
| 484 | * \retval true if the vertex is locked | ||
| 485 | * \retval false otherwise | ||
| 486 | */ | ||
| 487 | ✗ | bool is_locked(index_t v) { | |
| 488 | ✗ | return (v==locked_1_ || v==locked_2_); | |
| 489 | } | ||
| 490 | |||
| 491 | /** | ||
| 492 | * \brief Copies u,v coordinates from the mesh to OpenNL solver. | ||
| 493 | */ | ||
| 494 | ✗ | void mesh_to_solver() { | |
| 495 | ✗ | for(NLuint i=0; i<mesh_.vertices.nb(); ++i) { | |
| 496 | ✗ | double u = tex_coord_[2*i]; | |
| 497 | ✗ | double v = tex_coord_[2*i+1]; | |
| 498 | ✗ | nlSetVariable(2 * i , u); | |
| 499 | ✗ | nlSetVariable(2 * i + 1, v); | |
| 500 | ✗ | if(!spectral_ && is_locked(i)) { | |
| 501 | ✗ | nlLockVariable(2 * i ); | |
| 502 | ✗ | nlLockVariable(2 * i + 1); | |
| 503 | } | ||
| 504 | } | ||
| 505 | ✗ | } | |
| 506 | |||
| 507 | /** | ||
| 508 | * \brief Chooses an initial solution, and locks two vertices. | ||
| 509 | */ | ||
| 510 | ✗ | void project() { | |
| 511 | // Get bbox | ||
| 512 | ✗ | double xmin = 1e30; | |
| 513 | ✗ | double ymin = 1e30; | |
| 514 | ✗ | double zmin = 1e30; | |
| 515 | ✗ | double xmax = -1e30; | |
| 516 | ✗ | double ymax = -1e30; | |
| 517 | ✗ | double zmax = -1e30; | |
| 518 | |||
| 519 | ✗ | for(index_t i: mesh_.vertices) { | |
| 520 | ✗ | const vec3& p = mesh_.vertices.point(i); | |
| 521 | ✗ | xmin = std::min(p.x, xmin); | |
| 522 | ✗ | ymin = std::min(p.y, ymin); | |
| 523 | ✗ | zmin = std::min(p.z, zmin); | |
| 524 | |||
| 525 | ✗ | xmax = std::max(p.x, xmax); | |
| 526 | ✗ | ymax = std::max(p.y, ymax); | |
| 527 | ✗ | zmax = std::max(p.z, zmax); | |
| 528 | } | ||
| 529 | |||
| 530 | ✗ | double dx = xmax - xmin; | |
| 531 | ✗ | double dy = ymax - ymin; | |
| 532 | ✗ | double dz = zmax - zmin; | |
| 533 | |||
| 534 | ✗ | vec3 V1,V2; | |
| 535 | |||
| 536 | // Find shortest bbox axis | ||
| 537 | ✗ | if(dx <= dy && dx <= dz) { | |
| 538 | ✗ | if(dy > dz) { | |
| 539 | ✗ | V1 = vec3(0,1,0); | |
| 540 | ✗ | V2 = vec3(0,0,1); | |
| 541 | } else { | ||
| 542 | ✗ | V2 = vec3(0,1,0); | |
| 543 | ✗ | V1 = vec3(0,0,1); | |
| 544 | } | ||
| 545 | ✗ | } else if(dy <= dx && dy <= dz) { | |
| 546 | ✗ | if(dx > dz) { | |
| 547 | ✗ | V1 = vec3(1,0,0); | |
| 548 | ✗ | V2 = vec3(0,0,1); | |
| 549 | } else { | ||
| 550 | ✗ | V2 = vec3(1,0,0); | |
| 551 | ✗ | V1 = vec3(0,0,1); | |
| 552 | } | ||
| 553 | ✗ | } else if(dz <= dx && dz <= dy) { | |
| 554 | ✗ | if(dx > dy) { | |
| 555 | ✗ | V1 = vec3(1,0,0); | |
| 556 | ✗ | V2 = vec3(0,1,0); | |
| 557 | } else { | ||
| 558 | ✗ | V2 = vec3(1,0,0); | |
| 559 | ✗ | V1 = vec3(0,1,0); | |
| 560 | } | ||
| 561 | } | ||
| 562 | |||
| 563 | // Project onto shortest bbox axis, | ||
| 564 | // and lock extrema vertices | ||
| 565 | |||
| 566 | ✗ | double umin = 1e30; | |
| 567 | ✗ | double umax = -1e30; | |
| 568 | |||
| 569 | ✗ | for(index_t i: mesh_.vertices) { | |
| 570 | ✗ | const vec3& p = mesh_.vertices.point(i); | |
| 571 | ✗ | double u = dot(p,V1); | |
| 572 | ✗ | double v = dot(p,V2); | |
| 573 | ✗ | tex_coord_[2*i] = u; | |
| 574 | ✗ | tex_coord_[2*i+1] = v; | |
| 575 | ✗ | if(u < umin) { | |
| 576 | ✗ | locked_1_ = i; | |
| 577 | ✗ | umin = u; | |
| 578 | } | ||
| 579 | ✗ | if(u > umax) { | |
| 580 | ✗ | locked_2_ = i; | |
| 581 | ✗ | umax = u; | |
| 582 | } | ||
| 583 | } | ||
| 584 | ✗ | } | |
| 585 | |||
| 586 | Mesh& mesh_; | ||
| 587 | |||
| 588 | Attribute<double>& tex_coord_; | ||
| 589 | |||
| 590 | Attribute<double>& angle_; | ||
| 591 | /** | ||
| 592 | * \brief true if spectral mode is used, | ||
| 593 | * false if locked least squares mode is used. | ||
| 594 | */ | ||
| 595 | bool spectral_; | ||
| 596 | |||
| 597 | /** | ||
| 598 | * \brief In spectral mode, the index of the first | ||
| 599 | * non-zero eigenvalue. | ||
| 600 | */ | ||
| 601 | NLuint eigen_; | ||
| 602 | |||
| 603 | /** | ||
| 604 | * \brief The indices of the two locked vertices. | ||
| 605 | */ | ||
| 606 | index_t locked_1_, locked_2_; | ||
| 607 | |||
| 608 | bool verbose_; | ||
| 609 | }; | ||
| 610 | |||
| 611 | } | ||
| 612 | |||
| 613 | namespace GEO { | ||
| 614 | |||
| 615 | ✗ | void mesh_compute_LSCM( | |
| 616 | Mesh& M, const std::string& attribute_name, bool spectral, | ||
| 617 | const std::string& angle_attribute_name, | ||
| 618 | bool verbose | ||
| 619 | ) { | ||
| 620 | ✗ | Attribute<double> tex_coord; | |
| 621 | ✗ | tex_coord.bind_if_is_defined(M.vertices.attributes(), attribute_name); | |
| 622 | ✗ | if(tex_coord.is_bound() && tex_coord.dimension() != 2) { | |
| 623 | ✗ | Logger::err("LSCM") << "Attribute " << attribute_name | |
| 624 | ✗ | << " already exists in mesh with dimension " | |
| 625 | ✗ | << tex_coord.dimension() | |
| 626 | ✗ | << " (expected 2)" | |
| 627 | ✗ | << std::endl; | |
| 628 | ✗ | return; | |
| 629 | } | ||
| 630 | ✗ | if(!tex_coord.is_bound()) { | |
| 631 | ✗ | tex_coord.create_vector_attribute( | |
| 632 | M.vertices.attributes(),attribute_name,2 | ||
| 633 | ); | ||
| 634 | } | ||
| 635 | ✗ | Attribute<double> angle; | |
| 636 | ✗ | angle.bind_if_is_defined( | |
| 637 | M.facet_corners.attributes(), angle_attribute_name | ||
| 638 | ); | ||
| 639 | ✗ | LSCM lscm(M,tex_coord,angle); | |
| 640 | ✗ | lscm.set_spectral(spectral); | |
| 641 | ✗ | lscm.set_verbose(verbose); | |
| 642 | ✗ | lscm.apply(); | |
| 643 | ✗ | } | |
| 644 | } | ||
| 645 |