diff --git a/CMakeLists.txt b/CMakeLists.txt index 70326c1a7fd..80d5854ca2c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -449,6 +449,9 @@ option(WITH_MOD_FLUID "Enable Mantaflow Fluid Simulation Framework" ON) option(WITH_MOD_REMESH "Enable Remesh Modifier" ON) option(WITH_MOD_OCEANSIM "Enable Ocean Modifier" ON) +# UV solvers. +option(WITH_UV_SLIM "SLIM UV unwrapping solver" ON) + # Image format support option(WITH_IMAGE_OPENEXR "Enable OpenEXR Support (http://www.openexr.com)" ON) option(WITH_IMAGE_OPENJPEG "Enable OpenJpeg Support (http://www.openjpeg.org)" ON) diff --git a/build_files/cmake/config/blender_full.cmake b/build_files/cmake/config/blender_full.cmake index 5d0ddc005b1..f17c0b35258 100644 --- a/build_files/cmake/config/blender_full.cmake +++ b/build_files/cmake/config/blender_full.cmake @@ -40,6 +40,7 @@ set(WITH_LZO ON CACHE BOOL "" FORCE) set(WITH_MOD_FLUID ON CACHE BOOL "" FORCE) set(WITH_MOD_OCEANSIM ON CACHE BOOL "" FORCE) set(WITH_MOD_REMESH ON CACHE BOOL "" FORCE) +set(WITH_UV_SLIM ON CACHE BOOL "" FORCE) set(WITH_NANOVDB ON CACHE BOOL "" FORCE) set(WITH_OPENAL ON CACHE BOOL "" FORCE) set(WITH_OPENCOLLADA ON CACHE BOOL "" FORCE) diff --git a/build_files/cmake/config/blender_lite.cmake b/build_files/cmake/config/blender_lite.cmake index 3f4f1a9ca81..1fdc0b6fb11 100644 --- a/build_files/cmake/config/blender_lite.cmake +++ b/build_files/cmake/config/blender_lite.cmake @@ -46,6 +46,7 @@ set(WITH_LZO OFF CACHE BOOL "" FORCE) set(WITH_MOD_FLUID OFF CACHE BOOL "" FORCE) set(WITH_MOD_OCEANSIM OFF CACHE BOOL "" FORCE) set(WITH_MOD_REMESH OFF CACHE BOOL "" FORCE) +set(WITH_UV_SLIM OFF CACHE BOOL "" FORCE) set(WITH_NANOVDB OFF CACHE BOOL "" FORCE) set(WITH_OPENAL OFF CACHE BOOL "" FORCE) set(WITH_OPENCOLLADA OFF CACHE BOOL "" FORCE) diff --git a/build_files/cmake/config/blender_release.cmake b/build_files/cmake/config/blender_release.cmake index 81e8a7c85ec..4af6a307497 100644 --- a/build_files/cmake/config/blender_release.cmake +++ b/build_files/cmake/config/blender_release.cmake @@ -44,6 +44,7 @@ set(WITH_LZO ON CACHE BOOL "" FORCE) set(WITH_MOD_FLUID ON CACHE BOOL "" FORCE) set(WITH_MOD_OCEANSIM ON CACHE BOOL "" FORCE) set(WITH_MOD_REMESH ON CACHE BOOL "" FORCE) +set(WITH_UV_SLIM ON CACHE BOOL "" FORCE) set(WITH_NANOVDB ON CACHE BOOL "" FORCE) set(WITH_OPENAL ON CACHE BOOL "" FORCE) set(WITH_OPENCOLLADA ON CACHE BOOL "" FORCE) diff --git a/doc/doxygen/doxygen.intern.h b/doc/doxygen/doxygen.intern.h index 3690d4a0440..d5ac6971ba9 100644 --- a/doc/doxygen/doxygen.intern.h +++ b/doc/doxygen/doxygen.intern.h @@ -61,6 +61,9 @@ /** \defgroup intern_sky_model Sky Model * \ingroup intern */ +/** \defgroup intern_slim SLIM Solver for UV Unwrapping + * \ingroup intern */ + /** \defgroup intern_utf_conv UTF-8/16 Conversion (utfconv) * \ingroup intern */ diff --git a/doc/license/MPL2-license.txt b/doc/license/MPL2-license.txt new file mode 100644 index 00000000000..a612ad9813b --- /dev/null +++ b/doc/license/MPL2-license.txt @@ -0,0 +1,373 @@ +Mozilla Public License Version 2.0 +================================== + +1. Definitions +-------------- + +1.1. "Contributor" + means each individual or legal entity that creates, contributes to + the creation of, or owns Covered Software. + +1.2. "Contributor Version" + means the combination of the Contributions of others (if any) used + by a Contributor and that particular Contributor's Contribution. + +1.3. "Contribution" + means Covered Software of a particular Contributor. + +1.4. "Covered Software" + means Source Code Form to which the initial Contributor has attached + the notice in Exhibit A, the Executable Form of such Source Code + Form, and Modifications of such Source Code Form, in each case + including portions thereof. + +1.5. "Incompatible With Secondary Licenses" + means + + (a) that the initial Contributor has attached the notice described + in Exhibit B to the Covered Software; or + + (b) that the Covered Software was made available under the terms of + version 1.1 or earlier of the License, but not also under the + terms of a Secondary License. + +1.6. "Executable Form" + means any form of the work other than Source Code Form. + +1.7. "Larger Work" + means a work that combines Covered Software with other material, in + a separate file or files, that is not Covered Software. + +1.8. "License" + means this document. + +1.9. "Licensable" + means having the right to grant, to the maximum extent possible, + whether at the time of the initial grant or subsequently, any and + all of the rights conveyed by this License. + +1.10. "Modifications" + means any of the following: + + (a) any file in Source Code Form that results from an addition to, + deletion from, or modification of the contents of Covered + Software; or + + (b) any new file in Source Code Form that contains any Covered + Software. + +1.11. "Patent Claims" of a Contributor + means any patent claim(s), including without limitation, method, + process, and apparatus claims, in any patent Licensable by such + Contributor that would be infringed, but for the grant of the + License, by the making, using, selling, offering for sale, having + made, import, or transfer of either its Contributions or its + Contributor Version. + +1.12. "Secondary License" + means either the GNU General Public License, Version 2.0, the GNU + Lesser General Public License, Version 2.1, the GNU Affero General + Public License, Version 3.0, or any later versions of those + licenses. + +1.13. "Source Code Form" + means the form of the work preferred for making modifications. + +1.14. "You" (or "Your") + means an individual or a legal entity exercising rights under this + License. For legal entities, "You" includes any entity that + controls, is controlled by, or is under common control with You. For + purposes of this definition, "control" means (a) the power, direct + or indirect, to cause the direction or management of such entity, + whether by contract or otherwise, or (b) ownership of more than + fifty percent (50%) of the outstanding shares or beneficial + ownership of such entity. + +2. License Grants and Conditions +-------------------------------- + +2.1. Grants + +Each Contributor hereby grants You a world-wide, royalty-free, +non-exclusive license: + +(a) under intellectual property rights (other than patent or trademark) + Licensable by such Contributor to use, reproduce, make available, + modify, display, perform, distribute, and otherwise exploit its + Contributions, either on an unmodified basis, with Modifications, or + as part of a Larger Work; and + +(b) under Patent Claims of such Contributor to make, use, sell, offer + for sale, have made, import, and otherwise transfer either its + Contributions or its Contributor Version. + +2.2. Effective Date + +The licenses granted in Section 2.1 with respect to any Contribution +become effective for each Contribution on the date the Contributor first +distributes such Contribution. + +2.3. Limitations on Grant Scope + +The licenses granted in this Section 2 are the only rights granted under +this License. No additional rights or licenses will be implied from the +distribution or licensing of Covered Software under this License. +Notwithstanding Section 2.1(b) above, no patent license is granted by a +Contributor: + +(a) for any code that a Contributor has removed from Covered Software; + or + +(b) for infringements caused by: (i) Your and any other third party's + modifications of Covered Software, or (ii) the combination of its + Contributions with other software (except as part of its Contributor + Version); or + +(c) under Patent Claims infringed by Covered Software in the absence of + its Contributions. + +This License does not grant any rights in the trademarks, service marks, +or logos of any Contributor (except as may be necessary to comply with +the notice requirements in Section 3.4). + +2.4. Subsequent Licenses + +No Contributor makes additional grants as a result of Your choice to +distribute the Covered Software under a subsequent version of this +License (see Section 10.2) or under the terms of a Secondary License (if +permitted under the terms of Section 3.3). + +2.5. Representation + +Each Contributor represents that the Contributor believes its +Contributions are its original creation(s) or it has sufficient rights +to grant the rights to its Contributions conveyed by this License. + +2.6. Fair Use + +This License is not intended to limit any rights You have under +applicable copyright doctrines of fair use, fair dealing, or other +equivalents. + +2.7. Conditions + +Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted +in Section 2.1. + +3. Responsibilities +------------------- + +3.1. Distribution of Source Form + +All distribution of Covered Software in Source Code Form, including any +Modifications that You create or to which You contribute, must be under +the terms of this License. You must inform recipients that the Source +Code Form of the Covered Software is governed by the terms of this +License, and how they can obtain a copy of this License. You may not +attempt to alter or restrict the recipients' rights in the Source Code +Form. + +3.2. Distribution of Executable Form + +If You distribute Covered Software in Executable Form then: + +(a) such Covered Software must also be made available in Source Code + Form, as described in Section 3.1, and You must inform recipients of + the Executable Form how they can obtain a copy of such Source Code + Form by reasonable means in a timely manner, at a charge no more + than the cost of distribution to the recipient; and + +(b) You may distribute such Executable Form under the terms of this + License, or sublicense it under different terms, provided that the + license for the Executable Form does not attempt to limit or alter + the recipients' rights in the Source Code Form under this License. + +3.3. Distribution of a Larger Work + +You may create and distribute a Larger Work under terms of Your choice, +provided that You also comply with the requirements of this License for +the Covered Software. If the Larger Work is a combination of Covered +Software with a work governed by one or more Secondary Licenses, and the +Covered Software is not Incompatible With Secondary Licenses, this +License permits You to additionally distribute such Covered Software +under the terms of such Secondary License(s), so that the recipient of +the Larger Work may, at their option, further distribute the Covered +Software under the terms of either this License or such Secondary +License(s). + +3.4. Notices + +You may not remove or alter the substance of any license notices +(including copyright notices, patent notices, disclaimers of warranty, +or limitations of liability) contained within the Source Code Form of +the Covered Software, except that You may alter any license notices to +the extent required to remedy known factual inaccuracies. + +3.5. Application of Additional Terms + +You may choose to offer, and to charge a fee for, warranty, support, +indemnity or liability obligations to one or more recipients of Covered +Software. However, You may do so only on Your own behalf, and not on +behalf of any Contributor. You must make it absolutely clear that any +such warranty, support, indemnity, or liability obligation is offered by +You alone, and You hereby agree to indemnify every Contributor for any +liability incurred by such Contributor as a result of warranty, support, +indemnity or liability terms You offer. You may include additional +disclaimers of warranty and limitations of liability specific to any +jurisdiction. + +4. Inability to Comply Due to Statute or Regulation +--------------------------------------------------- + +If it is impossible for You to comply with any of the terms of this +License with respect to some or all of the Covered Software due to +statute, judicial order, or regulation then You must: (a) comply with +the terms of this License to the maximum extent possible; and (b) +describe the limitations and the code they affect. Such description must +be placed in a text file included with all distributions of the Covered +Software under this License. Except to the extent prohibited by statute +or regulation, such description must be sufficiently detailed for a +recipient of ordinary skill to be able to understand it. + +5. Termination +-------------- + +5.1. The rights granted under this License will terminate automatically +if You fail to comply with any of its terms. However, if You become +compliant, then the rights granted under this License from a particular +Contributor are reinstated (a) provisionally, unless and until such +Contributor explicitly and finally terminates Your grants, and (b) on an +ongoing basis, if such Contributor fails to notify You of the +non-compliance by some reasonable means prior to 60 days after You have +come back into compliance. Moreover, Your grants from a particular +Contributor are reinstated on an ongoing basis if such Contributor +notifies You of the non-compliance by some reasonable means, this is the +first time You have received notice of non-compliance with this License +from such Contributor, and You become compliant prior to 30 days after +Your receipt of the notice. + +5.2. If You initiate litigation against any entity by asserting a patent +infringement claim (excluding declaratory judgment actions, +counter-claims, and cross-claims) alleging that a Contributor Version +directly or indirectly infringes any patent, then the rights granted to +You by any and all Contributors for the Covered Software under Section +2.1 of this License shall terminate. + +5.3. In the event of termination under Sections 5.1 or 5.2 above, all +end user license agreements (excluding distributors and resellers) which +have been validly granted by You or Your distributors under this License +prior to termination shall survive termination. + +************************************************************************ +* * +* 6. Disclaimer of Warranty * +* ------------------------- * +* * +* Covered Software is provided under this License on an "as is" * +* basis, without warranty of any kind, either expressed, implied, or * +* statutory, including, without limitation, warranties that the * +* Covered Software is free of defects, merchantable, fit for a * +* particular purpose or non-infringing. The entire risk as to the * +* quality and performance of the Covered Software is with You. * +* Should any Covered Software prove defective in any respect, You * +* (not any Contributor) assume the cost of any necessary servicing, * +* repair, or correction. This disclaimer of warranty constitutes an * +* essential part of this License. No use of any Covered Software is * +* authorized under this License except under this disclaimer. * +* * +************************************************************************ + +************************************************************************ +* * +* 7. Limitation of Liability * +* -------------------------- * +* * +* Under no circumstances and under no legal theory, whether tort * +* (including negligence), contract, or otherwise, shall any * +* Contributor, or anyone who distributes Covered Software as * +* permitted above, be liable to You for any direct, indirect, * +* special, incidental, or consequential damages of any character * +* including, without limitation, damages for lost profits, loss of * +* goodwill, work stoppage, computer failure or malfunction, or any * +* and all other commercial damages or losses, even if such party * +* shall have been informed of the possibility of such damages. This * +* limitation of liability shall not apply to liability for death or * +* personal injury resulting from such party's negligence to the * +* extent applicable law prohibits such limitation. Some * +* jurisdictions do not allow the exclusion or limitation of * +* incidental or consequential damages, so this exclusion and * +* limitation may not apply to You. * +* * +************************************************************************ + +8. Litigation +------------- + +Any litigation relating to this License may be brought only in the +courts of a jurisdiction where the defendant maintains its principal +place of business and such litigation shall be governed by laws of that +jurisdiction, without reference to its conflict-of-law provisions. +Nothing in this Section shall prevent a party's ability to bring +cross-claims or counter-claims. + +9. Miscellaneous +---------------- + +This License represents the complete agreement concerning the subject +matter hereof. If any provision of this License is held to be +unenforceable, such provision shall be reformed only to the extent +necessary to make it enforceable. Any law or regulation which provides +that the language of a contract shall be construed against the drafter +shall not be used to construe this License against a Contributor. + +10. Versions of the License +--------------------------- + +10.1. New Versions + +Mozilla Foundation is the license steward. Except as provided in Section +10.3, no one other than the license steward has the right to modify or +publish new versions of this License. Each version will be given a +distinguishing version number. + +10.2. Effect of New Versions + +You may distribute the Covered Software under the terms of the version +of the License under which You originally received the Covered Software, +or under the terms of any subsequent version published by the license +steward. + +10.3. Modified Versions + +If you create software not governed by this License, and you want to +create a new license for such software, you may create and use a +modified version of this License if you rename the license and remove +any references to the name of the license steward (except to note that +such modified license differs from this License). + +10.4. Distributing Source Code Form that is Incompatible With Secondary +Licenses + +If You choose to distribute Source Code Form that is Incompatible With +Secondary Licenses under the terms of this version of the License, the +notice described in Exhibit B of this License must be attached. + +Exhibit A - Source Code Form License Notice +------------------------------------------- + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + +If it is not possible or desirable to put the notice in a particular +file, then You may include the notice in a location (such as a LICENSE +file in a relevant directory) where a recipient would be likely to look +for such a notice. + +You may add additional accurate notices of copyright ownership. + +Exhibit B - "Incompatible With Secondary Licenses" Notice +--------------------------------------------------------- + + This Source Code Form is "Incompatible With Secondary Licenses", as + defined by the Mozilla Public License, v. 2.0. diff --git a/intern/CMakeLists.txt b/intern/CMakeLists.txt index 7b900d3beea..851e46077cf 100644 --- a/intern/CMakeLists.txt +++ b/intern/CMakeLists.txt @@ -52,6 +52,10 @@ if(WITH_MOD_FLUID) add_subdirectory(mantaflow) endif() +if(WITH_UV_SLIM) + add_subdirectory(slim) +endif() + if(WITH_OPENVDB) add_subdirectory(openvdb) endif() diff --git a/intern/slim/CMakeLists.txt b/intern/slim/CMakeLists.txt new file mode 100644 index 00000000000..55ac25447a3 --- /dev/null +++ b/intern/slim/CMakeLists.txt @@ -0,0 +1,42 @@ +# SPDX-FileCopyrightText: 2023 Blender Authors +# +# SPDX-License-Identifier: GPL-2.0-or-later + +# Derived from `libigl`, a simple C++ geometry processing library. + +set(INC + ./ + ../../source/blender/blenlib + ../../intern/guardedalloc +) + +set(INC_SYS + ${EIGEN3_INCLUDE_DIRS} +) + +set(SRC + slim_matrix_transfer.h + + intern/area_compensation.cpp + intern/area_compensation.h + intern/cotmatrix.h + intern/doublearea.h + intern/edge_lengths.h + intern/flip_avoiding_line_search.h + intern/geometry_data_retrieval.cpp + intern/geometry_data_retrieval.h + intern/least_squares_relocator.cpp + intern/least_squares_relocator.h + intern/slim.cpp + intern/slim.h + intern/slim_matrix_transfer.cpp + intern/slim_parametrizer.cpp + intern/uv_initializer.cpp + intern/uv_initializer.h +) + +set(LIB +) + +blender_add_lib(bf_intern_slim "${SRC}" "${INC}" "${INC_SYS}" "${LIB}") +add_library(bf::intern::slim ALIAS bf_intern_slim) diff --git a/intern/slim/intern/area_compensation.cpp b/intern/slim/intern/area_compensation.cpp new file mode 100644 index 00000000000..ff75affad81 --- /dev/null +++ b/intern/slim/intern/area_compensation.cpp @@ -0,0 +1,94 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#include "BLI_assert.h" + +#include + +#include "area_compensation.h" +#include "doublearea.h" +#include "slim.h" + +using namespace Eigen; + +namespace slim { + +static void correct_geometry_size(double surface_area_to_map_area_ratio, + MatrixXd &vertex_positions, + double desired_surface_area_to_map_ration) +{ + BLI_assert(surface_area_to_map_area_ratio > 0); + double sqrt_of_ratio = sqrt(surface_area_to_map_area_ratio / desired_surface_area_to_map_ration); + vertex_positions = vertex_positions / sqrt_of_ratio; +} + +template +static double compute_surface_area(const VertexPositionType v, const FaceIndicesType f) +{ + Eigen::VectorXd doubled_area_of_triangles; + doublearea(v, f, doubled_area_of_triangles); + double area_of_map = doubled_area_of_triangles.sum() / 2; + return area_of_map; +} + +void correct_map_surface_area_if_necessary(SLIMData &slim_data) +{ + if (!slim_data.valid) { + return; + } + + bool mesh_surface_area_was_corrected = (slim_data.expectedSurfaceAreaOfResultingMap != 0); + int number_of_pinned_vertices = slim_data.b.rows(); + bool no_pinned_vertices_exist = number_of_pinned_vertices == 0; + + bool needs_area_correction = mesh_surface_area_was_corrected && no_pinned_vertices_exist; + if (!needs_area_correction) { + return; + } + + double area_ofresulting_map = compute_surface_area(slim_data.V_o, slim_data.F); + if (!area_ofresulting_map) { + return; + } + + double resulting_area_to_expected_area_ratio = area_ofresulting_map / + slim_data.expectedSurfaceAreaOfResultingMap; + double desired_ratio = 1.0; + correct_geometry_size(resulting_area_to_expected_area_ratio, slim_data.V_o, desired_ratio); +} + +void correct_mesh_surface_area_if_necessary(SLIMData &slim_data) +{ + BLI_assert(slim_data.valid); + + int number_of_pinned_vertices = slim_data.b.rows(); + bool pinned_vertices_exist = number_of_pinned_vertices > 0; + bool needs_area_correction = slim_data.skipInitialization || pinned_vertices_exist; + + if (!needs_area_correction) { + return; + } + + double area_of_preinitialized_map = compute_surface_area(slim_data.V_o, slim_data.F); + if (!area_of_preinitialized_map) { + return; + } + + if (area_of_preinitialized_map < 0) { + area_of_preinitialized_map *= -1; + } + + slim_data.expectedSurfaceAreaOfResultingMap = area_of_preinitialized_map; + double surface_area_of3d_mesh = compute_surface_area(slim_data.V, slim_data.F); + double surface_area_to_map_area_ratio = surface_area_of3d_mesh / area_of_preinitialized_map; + + double desired_ratio = 1.0; + correct_geometry_size(surface_area_to_map_area_ratio, slim_data.V, desired_ratio); +} + +} // namespace slim diff --git a/intern/slim/intern/area_compensation.h b/intern/slim/intern/area_compensation.h new file mode 100644 index 00000000000..28b42b3ab02 --- /dev/null +++ b/intern/slim/intern/area_compensation.h @@ -0,0 +1,22 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include "slim.h" + +#include + +using namespace Eigen; + +namespace slim { + +void correct_map_surface_area_if_necessary(SLIMData &slimData); +void correct_mesh_surface_area_if_necessary(SLIMData &slimData); + +} // namespace slim diff --git a/intern/slim/intern/cotmatrix.cpp b/intern/slim/intern/cotmatrix.cpp new file mode 100644 index 00000000000..d7c0d8cb5fe --- /dev/null +++ b/intern/slim/intern/cotmatrix.cpp @@ -0,0 +1,97 @@ +/* SPDX-FileCopyrightText: 2013 Alec Jacobson + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#include "cotmatrix.h" +#include "edge_lengths.h" + +#include + +namespace slim { + +/* Inputs: + * V #V by dim list of rest domain positions + * F #F by 3 list of triangle indices into V + * Outputs: + * C #F by 3 list of 1/2*cotangents corresponding angles + * for triangles, columns correspond to edges [1,2],[2,0],[0,1] + */ +template +static inline void cotmatrix_entries(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &C) +{ + using namespace std; + using namespace Eigen; + /* Number of elements. */ + int m = F.rows(); + + assert(F.cols() == 3); + + /* Law of cosines + law of sines. */ + /* Compute Squared Edge lenghts. */ + Matrix l2; + squared_edge_lengths(V, F, l2); + /* Compute Edge lenghts. */ + Matrix l; + l = l2.array().sqrt(); + + /* Double area. */ + Matrix dblA; + doublearea(l, dblA); + /* Cotangents and diagonal entries for element matrices. + * correctly divided by 4. */ + C.resize(m, 3); + for (int i = 0; i < m; i++) { + C(i, 0) = (l2(i, 1) + l2(i, 2) - l2(i, 0)) / dblA(i) / 4.0; + C(i, 1) = (l2(i, 2) + l2(i, 0) - l2(i, 1)) / dblA(i) / 4.0; + C(i, 2) = (l2(i, 0) + l2(i, 1) - l2(i, 2)) / dblA(i) / 4.0; + } +} + +template +inline void cotmatrix(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::SparseMatrix &L) +{ + using namespace Eigen; + using namespace std; + + L.resize(V.rows(), V.rows()); + Matrix edges; + /* 3 for triangles. */ + assert(F.cols() == 3); + + /* This is important! it could decrease the comptuation time by a factor of 2 + * Laplacian for a closed 2d manifold mesh will have on average 7 entries per + * row. */ + L.reserve(10 * V.rows()); + edges.resize(3, 2); + edges << 1, 2, 2, 0, 0, 1; + + /* Gather cotangents. */ + Matrix C; + cotmatrix_entries(V, F, C); + + vector> IJV; + IJV.reserve(F.rows() * edges.rows() * 4); + /* Loop over triangles. */ + for (int i = 0; i < F.rows(); i++) { + /* Loop over edges of element. */ + for (int e = 0; e < edges.rows(); e++) { + int source = F(i, edges(e, 0)); + int dest = F(i, edges(e, 1)); + IJV.push_back(Triplet(source, dest, C(i, e))); + IJV.push_back(Triplet(dest, source, C(i, e))); + IJV.push_back(Triplet(source, source, -C(i, e))); + IJV.push_back(Triplet(dest, dest, -C(i, e))); + } + } + L.setFromTriplets(IJV.begin(), IJV.end()); +} + +} // namespace slim diff --git a/intern/slim/intern/cotmatrix.h b/intern/slim/intern/cotmatrix.h new file mode 100644 index 00000000000..1ce5fafdada --- /dev/null +++ b/intern/slim/intern/cotmatrix.h @@ -0,0 +1,44 @@ +/* SPDX-FileCopyrightText: 2014 Alec Jacobson + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include +#include + +namespace slim { + +/* Constructs the cotangent stiffness matrix (discrete laplacian) for a given + * mesh (V,F). + * + * Templates: + * DerivedV derived type of eigen matrix for V (e.g. derived from + * MatrixXd) + * DerivedF derived type of eigen matrix for F (e.g. derived from + * MatrixXi) + * Scalar scalar type for eigen sparse matrix (e.g. double) + * Inputs: + * V #V by dim list of mesh vertex positions + * F #F by simplex_size list of mesh faces (must be triangles) + * Outputs: + * L #V by #V cotangent matrix, each row i corresponding to V(i,:) + * + * Note: This Laplacian uses the convention that diagonal entries are + * **minus** the sum of off-diagonal entries. The diagonal entries are + * therefore in general negative and the matrix is **negative** semi-definite + * (immediately, -L is **positive** semi-definite) + */ +template +inline void cotmatrix(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::SparseMatrix &L); + +} // namespace slim + +#include "cotmatrix.cpp" diff --git a/intern/slim/intern/doublearea.cpp b/intern/slim/intern/doublearea.cpp new file mode 100644 index 00000000000..c8af534d6f0 --- /dev/null +++ b/intern/slim/intern/doublearea.cpp @@ -0,0 +1,199 @@ +/* SPDX-FileCopyrightText: 2013 Alec Jacobson + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#include "doublearea.h" +#include "edge_lengths.h" + +#include + +#include + +namespace slim { + +/* Sort the elements of a matrix X along a given dimension like matlabs sort + * function, assuming X.cols() == 3. + * + * Templates: + * DerivedX derived scalar type, e.g. MatrixXi or MatrixXd + * DerivedIX derived integer type, e.g. MatrixXi + * Inputs: + * X m by n matrix whose entries are to be sorted + * dim dimensional along which to sort: + * 1 sort each column (matlab default) + * 2 sort each row + * ascending sort ascending (true, matlab default) or descending (false) + * Outputs: + * Y m by n matrix whose entries are sorted + * IX m by n matrix of indices so that if dim = 1, then in matlab notation + * for j = 1:n, Y(:,j) = X(I(:,j),j); end + */ +template +static inline void doublearea_sort3(const Eigen::PlainObjectBase &X, + const int dim, + const bool ascending, + Eigen::PlainObjectBase &Y, + Eigen::PlainObjectBase &IX) +{ + using namespace Eigen; + using namespace std; + typedef typename Eigen::PlainObjectBase::Scalar YScalar; + Y = X.template cast(); + /* Get number of columns (or rows). */ + int num_outer = (dim == 1 ? X.cols() : X.rows()); + /* Get number of rows (or columns). */ + int num_inner = (dim == 1 ? X.rows() : X.cols()); + assert(num_inner == 3); + (void)num_inner; + typedef typename Eigen::PlainObjectBase::Scalar Index; + IX.resize(X.rows(), X.cols()); + if (dim == 1) { + IX.row(0).setConstant(0); /* = Eigen::PlainObjectBase::Zero(1,IX.cols());. */ + IX.row(1).setConstant(1); /* = Eigen::PlainObjectBase::Ones (1,IX.cols());. */ + IX.row(2).setConstant(2); /* = Eigen::PlainObjectBase::Ones (1,IX.cols());. */ + } + else { + IX.col(0).setConstant(0); /* = Eigen::PlainObjectBase::Zero(IX.rows(),1);. */ + IX.col(1).setConstant(1); /* = Eigen::PlainObjectBase::Ones (IX.rows(),1);. */ + IX.col(2).setConstant(2); /* = Eigen::PlainObjectBase::Ones (IX.rows(),1);. */ + } + + using namespace blender; + threading::parallel_for( + IndexRange(num_outer), 16000, [&IX, &Y, &dim, &ascending](const IndexRange range) { + for (const Index i : range) { + YScalar &a = (dim == 1 ? Y(0, i) : Y(i, 0)); + YScalar &b = (dim == 1 ? Y(1, i) : Y(i, 1)); + YScalar &c = (dim == 1 ? Y(2, i) : Y(i, 2)); + Index &ai = (dim == 1 ? IX(0, i) : IX(i, 0)); + Index &bi = (dim == 1 ? IX(1, i) : IX(i, 1)); + Index &ci = (dim == 1 ? IX(2, i) : IX(i, 2)); + if (ascending) { + /* 123 132 213 231 312 321. */ + if (a > b) { + std::swap(a, b); + std::swap(ai, bi); + } + /* 123 132 123 231 132 231. */ + if (b > c) { + std::swap(b, c); + std::swap(bi, ci); + /* 123 123 123 213 123 213. */ + if (a > b) { + std::swap(a, b); + std::swap(ai, bi); + } + /* 123 123 123 123 123 123. */ + } + } + else { + /* 123 132 213 231 312 321. */ + if (a < b) { + std::swap(a, b); + std::swap(ai, bi); + } + /* 213 312 213 321 312 321. */ + if (b < c) { + std::swap(b, c); + std::swap(bi, ci); + /* 231 321 231 321 321 321. */ + if (a < b) { + std::swap(a, b); + std::swap(ai, bi); + } + /* 321 321 321 321 321 321. */ + } + } + } + }); +} + +template +inline void doublearea(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &dblA) +{ + const int dim = V.cols(); + /* Only support triangles. */ + assert(F.cols() == 3); + const size_t m = F.rows(); + /* Compute edge lengths. */ + Eigen::Matrix l; + + /* Projected area helper. */ + const auto &proj_doublearea = [&V, &F](const int x, const int y, const int f) -> double { + auto rx = V(F(f, 0), x) - V(F(f, 2), x); + auto sx = V(F(f, 1), x) - V(F(f, 2), x); + auto ry = V(F(f, 0), y) - V(F(f, 2), y); + auto sy = V(F(f, 1), y) - V(F(f, 2), y); + return rx * sy - ry * sx; + }; + + switch (dim) { + case 3: { + dblA = Eigen::PlainObjectBase::Zero(m, 1); + for (size_t f = 0; f < m; f++) { + for (int d = 0; d < 3; d++) { + double dblAd = proj_doublearea(d, (d + 1) % 3, f); + dblA(f) += dblAd * dblAd; + } + } + dblA = dblA.array().sqrt().eval(); + break; + } + case 2: { + dblA.resize(m, 1); + for (size_t f = 0; f < m; f++) { + dblA(f) = proj_doublearea(0, 1, f); + } + break; + } + default: { + edge_lengths(V, F, l); + return doublearea(l, dblA); + } + } +} + +template +inline void doublearea(const Eigen::PlainObjectBase &ul, + Eigen::PlainObjectBase &dblA) +{ + using namespace Eigen; + using namespace std; + typedef typename Derivedl::Index Index; + /* Only support triangles. */ + assert(ul.cols() == 3); + /* Number of triangles. */ + const Index m = ul.rows(); + Eigen::Matrix l; + MatrixXi _; + /* "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1 + * http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf + * + * "Miscalculating Area and Angles of a Needle-like Triangle" + * https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf + */ + doublearea_sort3(ul, 2, false, l, _); + dblA.resize(l.rows(), 1); + + using namespace blender; + threading::parallel_for(IndexRange(m), 1000, [&l, &dblA](const IndexRange range) { + for (const Index i : range) { + /* Kahan's Heron's formula. */ + const typename Derivedl::Scalar arg = (l(i, 0) + (l(i, 1) + l(i, 2))) * + (l(i, 2) - (l(i, 0) - l(i, 1))) * + (l(i, 2) + (l(i, 0) - l(i, 1))) * + (l(i, 0) + (l(i, 1) - l(i, 2))); + dblA(i) = 2.0 * 0.25 * sqrt(arg); + assert(l(i, 2) - (l(i, 0) - l(i, 1)) && "FAILED KAHAN'S ASSERTION"); + assert(dblA(i) == dblA(i) && "DOUBLEAREA() PRODUCED NaN"); + } + }); +} + +} // namespace slim diff --git a/intern/slim/intern/doublearea.h b/intern/slim/intern/doublearea.h new file mode 100644 index 00000000000..dedb2edc5e1 --- /dev/null +++ b/intern/slim/intern/doublearea.h @@ -0,0 +1,52 @@ +/* SPDX-FileCopyrightText: 2013 Alec Jacobson + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include + +namespace slim { + +/* DOUBLEAREA computes twice the area for each input triangle + * + * Templates: + * DerivedV derived type of eigen matrix for V (e.g. derived from + * MatrixXd) + * DerivedF derived type of eigen matrix for F (e.g. derived from + * MatrixXi) + * DeriveddblA derived type of eigen matrix for dblA (e.g. derived from + * MatrixXd) + * Inputs: + * V #V by dim list of mesh vertex positions + * F #F by simplex_size list of mesh faces (must be triangles) + * Outputs: + * dblA #F list of triangle double areas (SIGNED only for 2D input) + * + * Known bug: For dim==3 complexity is O(#V + #F)!! Not just O(#F). This is a big deal + * if you have 1million unreferenced vertices and 1 face + */ +template +inline void doublearea(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &dblA); + +/* Same as above but use instrinsic edge lengths rather than (V,F) mesh + * Inputs: + * l #F by dim list of edge lengths using + * for triangles, columns correspond to edges 23,31,12 + * Outputs: + * dblA #F list of triangle double areas + */ +template +inline void doublearea(const Eigen::PlainObjectBase &l, + Eigen::PlainObjectBase &dblA); + +} // namespace slim + +#include "doublearea.cpp" diff --git a/intern/slim/intern/edge_lengths.cpp b/intern/slim/intern/edge_lengths.cpp new file mode 100644 index 00000000000..95f0340185b --- /dev/null +++ b/intern/slim/intern/edge_lengths.cpp @@ -0,0 +1,47 @@ +/* SPDX-FileCopyrightText: 2013 Alec Jacobson + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#include "BLI_task.hh" + +#include "edge_lengths.h" + +namespace slim { + +template +inline void edge_lengths(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &L) +{ + squared_edge_lengths(V, F, L); + L = L.array().sqrt().eval(); +} + +template +inline void squared_edge_lengths(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &L) +{ + using namespace std; + const int m = F.rows(); + assert(F.cols() == 3); + + L.resize(m, 3); + + /* Loop over faces. */ + using namespace blender; + threading::parallel_for(IndexRange(m), 1000, [&V, &F, &L](const IndexRange range) { + for (const int i : range) { + L(i, 0) = (V.row(F(i, 1)) - V.row(F(i, 2))).squaredNorm(); + L(i, 1) = (V.row(F(i, 2)) - V.row(F(i, 0))).squaredNorm(); + L(i, 2) = (V.row(F(i, 0)) - V.row(F(i, 1))).squaredNorm(); + } + }); +} + +} // namespace slim diff --git a/intern/slim/intern/edge_lengths.h b/intern/slim/intern/edge_lengths.h new file mode 100644 index 00000000000..a68c003c944 --- /dev/null +++ b/intern/slim/intern/edge_lengths.h @@ -0,0 +1,62 @@ +/* SPDX-FileCopyrightText: 2013 Alec Jacobson + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include + +namespace slim { +/* Constructs a list of lengths of edges opposite each index in a face + * (triangle) list + * + * Templates: + * DerivedV derived from vertex positions matrix type: i.e. MatrixXd + * DerivedF derived from face indices matrix type: i.e. MatrixXi + * DerivedL derived from edge lengths matrix type: i.e. MatrixXd + * Inputs: + * V eigen matrix #V by 3 + * F #F by 2 list of mesh edges + * or + * F #F by 3 list of mesh faces (must be triangles) + * or + * T #T by 4 list of mesh elements (must be tets) + * Outputs: + * L #F by {1|3|6} list of edge lengths + * for edges, column of lengths + * for triangles, columns correspond to edges [1,2],[2,0],[0,1] + */ +template +inline void edge_lengths(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &L); + +/* Constructs a list of squared lengths of edges opposite each index in a face + * (triangle) list + * + * Templates: + * DerivedV derived from vertex positions matrix type: i.e. MatrixXd + * DerivedF derived from face indices matrix type: i.e. MatrixXi + * DerivedL derived from edge lengths matrix type: i.e. MatrixXd + * Inputs: + * V eigen matrix #V by 3 + * F #F by 2 list of mesh edges + * or + * F #F by 3 list of mesh faces (must be triangles) + * Outputs: + * L #F by {1|3|6} list of edge lengths squared + * for edges, column of lengths + * for triangles, columns correspond to edges [1,2],[2,0],[0,1] + */ +template +inline void squared_edge_lengths(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &L); +} // namespace slim + +#include "edge_lengths.cpp" diff --git a/intern/slim/intern/flip_avoiding_line_search.cpp b/intern/slim/intern/flip_avoiding_line_search.cpp new file mode 100644 index 00000000000..617fb512384 --- /dev/null +++ b/intern/slim/intern/flip_avoiding_line_search.cpp @@ -0,0 +1,169 @@ +/* SPDX-FileCopyrightText: 2016 Michael Rabinovich + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#include "flip_avoiding_line_search.h" + +#include +#include + +namespace slim { + +/* Implement a bisection linesearch to minimize a mesh-based energy on vertices given at 'x' at a + * search direction 'd', with initial step size. Stops when a point with lower energy is found, or + * after maximal iterations have been reached. + * + * Inputs: + * x #X by dim list of variables + * d #X by dim list of a given search direction + * step_size initial step size + * energy A function to compute the mesh-based energy (return an energy that is + * bigger than 0) cur_energy(OPTIONAL) The energy at the given point. Helps save redundant + * computations. + * This is optional. If not specified, the function will compute it. + * Outputs: + * x #X by dim list of variables at the new location + * Returns the energy at the new point 'x'. + */ +static inline double line_search(Eigen::MatrixXd &x, + const Eigen::MatrixXd &d, + double step_size, + std::function energy, + double cur_energy = -1) +{ + double old_energy; + if (cur_energy > 0) { + old_energy = cur_energy; + } + else { + old_energy = energy(x); /* No energy was given -> need to compute the current energy. */ + } + double new_energy = old_energy; + int cur_iter = 0; + int MAX_STEP_SIZE_ITER = 12; + + while (new_energy >= old_energy && cur_iter < MAX_STEP_SIZE_ITER) { + Eigen::MatrixXd new_x = x + step_size * d; + + double cur_e = energy(new_x); + if (cur_e >= old_energy) { + step_size /= 2; + } + else { + x = new_x; + new_energy = cur_e; + } + cur_iter++; + } + return new_energy; +} + +static inline double get_smallest_pos_quad_zero(double a, double b, double c) +{ + using namespace std; + double t1, t2; + if (a != 0) { + double delta_in = pow(b, 2) - 4 * a * c; + if (delta_in < 0) { + return INFINITY; + } + double delta = sqrt(delta_in); + t1 = (-b + delta) / (2 * a); + t2 = (-b - delta) / (2 * a); + } + else { + t1 = t2 = -b / c; + } + + if (!std::isfinite(t1) || !std::isfinite(t2)) { + throw SlimFailedException(); + } + + double tmp_n = min(t1, t2); + t1 = max(t1, t2); + t2 = tmp_n; + if (t1 == t2) { + return INFINITY; /* Means the orientation flips twice = doesn't flip. */ + } + /* Return the smallest negative root if it exists, otherwise return infinity. */ + if (t1 > 0) { + if (t2 > 0) { + return t2; + } + else { + return t1; + } + } + else { + return INFINITY; + } +} + +static inline double get_min_pos_root_2D(const Eigen::MatrixXd &uv, + const Eigen::MatrixXi &F, + Eigen::MatrixXd &d, + int f) +{ + using namespace std; + /* Finding the smallest timestep t s.t a triangle get degenerated (<=> det = 0). */ + int v1 = F(f, 0); + int v2 = F(f, 1); + int v3 = F(f, 2); + /* Get quadratic coefficients (ax^2 + b^x + c). */ + const double &U11 = uv(v1, 0); + const double &U12 = uv(v1, 1); + const double &U21 = uv(v2, 0); + const double &U22 = uv(v2, 1); + const double &U31 = uv(v3, 0); + const double &U32 = uv(v3, 1); + + const double &V11 = d(v1, 0); + const double &V12 = d(v1, 1); + const double &V21 = d(v2, 0); + const double &V22 = d(v2, 1); + const double &V31 = d(v3, 0); + const double &V32 = d(v3, 1); + + double a = V11 * V22 - V12 * V21 - V11 * V32 + V12 * V31 + V21 * V32 - V22 * V31; + double b = U11 * V22 - U12 * V21 - U21 * V12 + U22 * V11 - U11 * V32 + U12 * V31 + U31 * V12 - + U32 * V11 + U21 * V32 - U22 * V31 - U31 * V22 + U32 * V21; + double c = U11 * U22 - U12 * U21 - U11 * U32 + U12 * U31 + U21 * U32 - U22 * U31; + + return get_smallest_pos_quad_zero(a, b, c); +} + +static inline double compute_max_step_from_singularities(const Eigen::MatrixXd &uv, + const Eigen::MatrixXi &F, + Eigen::MatrixXd &d) +{ + using namespace std; + double max_step = INFINITY; + + /* The if statement is outside the for loops to avoid branching/ease parallelizing. */ + for (int f = 0; f < F.rows(); f++) { + double min_positive_root = get_min_pos_root_2D(uv, F, d, f); + max_step = min(max_step, min_positive_root); + } + return max_step; +} + +inline double flip_avoiding_line_search(const Eigen::MatrixXi F, + Eigen::MatrixXd &cur_v, + Eigen::MatrixXd &dst_v, + std::function energy, + double cur_energy) +{ + using namespace std; + + Eigen::MatrixXd d = dst_v - cur_v; + double min_step_to_singularity = compute_max_step_from_singularities(cur_v, F, d); + double max_step_size = min(1., min_step_to_singularity * 0.8); + return line_search(cur_v, d, max_step_size, energy, cur_energy); +} + +} // namespace slim diff --git a/intern/slim/intern/flip_avoiding_line_search.h b/intern/slim/intern/flip_avoiding_line_search.h new file mode 100644 index 00000000000..0b24cc79f54 --- /dev/null +++ b/intern/slim/intern/flip_avoiding_line_search.h @@ -0,0 +1,46 @@ +/* SPDX-FileCopyrightText: 2016 Michael Rabinovich + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include + +namespace slim { + +/* A bisection line search for a mesh based energy that avoids triangle flips as suggested in + * "Bijective Parameterization with Free Boundaries" (Smith J. and Schaefer S., 2015). + * + * The user specifies an initial vertices position (that has no flips) and target one (that my have + * flipped triangles). This method first computes the largest step in direction of the destination + * vertices that does not incur flips, and then minimizes a given energy using this maximal step + * and a bisection linesearch (see igl::line_search). + * + * Supports triangle meshes. + * + * Inputs: + * F #F by 3 list of mesh faces + * cur_v #V by dim list of variables + * dst_v #V by dim list of target vertices. This mesh may have flipped triangles + * energy A function to compute the mesh-based energy (return an energy that is + * bigger than 0) cur_energy(OPTIONAL) The energy at the given point. Helps save + * redundant computations. + * This is optional. If not specified, the function will compute it. + * Outputs: + * cur_v #V by dim list of variables at the new location + * Returns the energy at the new point. + */ +inline double flip_avoiding_line_search(const Eigen::MatrixXi F, + Eigen::MatrixXd &cur_v, + Eigen::MatrixXd &dst_v, + std::function energy, + double cur_energy = -1); + +} // namespace slim + +#include "flip_avoiding_line_search.cpp" diff --git a/intern/slim/intern/geometry_data_retrieval.cpp b/intern/slim/intern/geometry_data_retrieval.cpp new file mode 100644 index 00000000000..8173caf41ed --- /dev/null +++ b/intern/slim/intern/geometry_data_retrieval.cpp @@ -0,0 +1,227 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#include "geometry_data_retrieval.h" + +#include "slim.h" +#include "slim_matrix_transfer.h" +#include + +#include "BLI_assert.h" +#include "area_compensation.h" +#include "geometry_data_retrieval.h" +#include "least_squares_relocator.h" + +#include "BLI_assert.h" + +using namespace Eigen; + +namespace slim { + +GeometryData::GeometryData(const MatrixTransfer &mt, MatrixTransferChart &chart) + : number_of_vertices(chart.verts_num), + number_of_faces(chart.faces_num), + /* `n_edges` in transferred_data accounts for boundary edges only once. */ + number_of_edges_twice(chart.edges_num + chart.boundary_vertices_num), + number_of_boundary_vertices(chart.boundary_vertices_num), + number_of_pinned_vertices(chart.pinned_vertices_num), + use_weights(mt.use_weights), + weight_influence(mt.weight_influence), + vertex_positions3d(chart.v_matrices.data(), number_of_vertices, columns_3), + uv_positions2d(chart.uv_matrices.data(), number_of_vertices, columns_2), + positions_of_pinned_vertices2d(), + positions_of_explicitly_pinned_vertices2d( + number_of_pinned_vertices != 0 ? chart.pp_matrices.data() : nullptr, + number_of_pinned_vertices, + columns_2), + faces_by_vertexindices(chart.f_matrices.data(), number_of_faces, columns_3), + edges_by_vertexindices(chart.e_matrices.data(), number_of_edges_twice, columns_2), + pinned_vertex_indices(), + explicitly_pinned_vertex_indices(number_of_pinned_vertices != 0 ? chart.p_matrices.data() : + nullptr, + number_of_pinned_vertices), + edge_lengths(chart.el_vectors.data(), number_of_edges_twice), + boundary_vertex_indices(chart.b_vectors.data(), number_of_boundary_vertices), + weights_per_vertex(chart.w_vectors.data(), number_of_vertices) +{ + retrieve_pinned_vertices(mt.fixed_boundary); +} + +static void create_weights_per_face(SLIMData &slim_data) +{ + if (!slim_data.valid) { + return; + } + + if (!slim_data.withWeightedParameterization) { + slim_data.weightPerFaceMap = Eigen::VectorXf::Ones(slim_data.F.rows()); + return; + } + + slim_data.weightPerFaceMap = Eigen::VectorXf(slim_data.F.rows()); + + /* The actual weight is `max_factor ^ (2 * (mean - 0.5))` */ + int weight_influence_sign = (slim_data.weightInfluence >= 0) ? 1 : -1; + double max_factor = std::abs(slim_data.weightInfluence) + 1; + + for (int fid = 0; fid < slim_data.F.rows(); fid++) { + Eigen::RowVector3i row = slim_data.F.row(fid); + float w1, w2, w3, mean, weight_factor, flipped_mean; + w1 = slim_data.weightmap(row(0)); + w2 = slim_data.weightmap(row(1)); + w3 = slim_data.weightmap(row(2)); + mean = (w1 + w2 + w3) / 3; + flipped_mean = 1 - mean; + + weight_factor = std::pow(max_factor, weight_influence_sign * 2 * (flipped_mean - 0.5)); + slim_data.weightPerFaceMap(fid) = weight_factor; + } +} + +void GeometryData::set_geometry_data_matrices(SLIMData &slim_data) const +{ + if (!slim_data.valid) { + return; + } + + slim_data.V = vertex_positions3d; + slim_data.F = faces_by_vertexindices; + slim_data.b = pinned_vertex_indices; + slim_data.bc = positions_of_pinned_vertices2d; + slim_data.V_o = uv_positions2d; + slim_data.oldUVs = uv_positions2d; + slim_data.weightmap = weights_per_vertex; + create_weights_per_face(slim_data); +} + +bool GeometryData::has_valid_preinitialized_map() const +{ + if (uv_positions2d.rows() == vertex_positions3d.rows() && uv_positions2d.cols() == columns_2) { + + int number_of_flips = count_flips(faces_by_vertexindices, uv_positions2d); + bool no_flips_present = (number_of_flips == 0); + return (no_flips_present); + } + return false; +} + +/* If we use interactive parametrisation, we usually start form an existing, flip-free unwrapping. + * Also, pinning of vertices has some issues with initialisation with convex border. + * We therefore may want to skip initialization. however, to skip initialization we need a + * preexisting valid starting map. */ +bool GeometryData::can_initialization_be_skipped(bool skip_initialization) const +{ + return (skip_initialization && has_valid_preinitialized_map()); +} + +void GeometryData::construct_slim_data(SLIMData &slim_data, + bool skip_initialization, + int reflection_mode) const +{ + BLI_assert(slim_data.valid); + + slim_data.skipInitialization = can_initialization_be_skipped(skip_initialization); + slim_data.weightInfluence = weight_influence; + slim_data.reflection_mode = reflection_mode; + slim_data.withWeightedParameterization = use_weights; + set_geometry_data_matrices(slim_data); + + double penalty_for_violating_pinned_positions = 10.0e100; + slim_data.soft_const_p = penalty_for_violating_pinned_positions; + slim_data.slim_energy = SLIMData::SYMMETRIC_DIRICHLET; + + initialize_if_needed(slim_data); + + transform_initialization_if_necessary(slim_data); + correct_mesh_surface_area_if_necessary(slim_data); + + slim_precompute(slim_data.V, + slim_data.F, + slim_data.V_o, + slim_data, + slim_data.slim_energy, + slim_data.b, + slim_data.bc, + slim_data.soft_const_p); +} + +void GeometryData::combine_matrices_of_pinned_and_boundary_vertices() +{ + /* Over - allocate pessimistically to avoid multiple reallocation. */ + int upper_bound_on_number_of_pinned_vertices = number_of_boundary_vertices + + number_of_pinned_vertices; + pinned_vertex_indices = VectorXi(upper_bound_on_number_of_pinned_vertices); + positions_of_pinned_vertices2d = MatrixXd(upper_bound_on_number_of_pinned_vertices, columns_2); + + /* Since border vertices use vertex indices 0 ... #bordervertices we can do: */ + pinned_vertex_indices.segment(0, number_of_boundary_vertices) = boundary_vertex_indices; + positions_of_pinned_vertices2d.block(0, 0, number_of_boundary_vertices, columns_2) = + uv_positions2d.block(0, 0, number_of_boundary_vertices, columns_2); + + int index = number_of_boundary_vertices; + int highest_vertex_index = (boundary_vertex_indices)(index - 1); + + for (Map::InnerIterator it(explicitly_pinned_vertex_indices, 0); it; ++it) { + int vertex_index = it.value(); + if (vertex_index > highest_vertex_index) { + pinned_vertex_indices(index) = vertex_index; + positions_of_pinned_vertices2d.row(index) = uv_positions2d.row(vertex_index); + index++; + } + } + + int actual_number_of_pinned_vertices = index; + pinned_vertex_indices.conservativeResize(actual_number_of_pinned_vertices); + positions_of_pinned_vertices2d.conservativeResize(actual_number_of_pinned_vertices, columns_2); + + number_of_pinned_vertices = actual_number_of_pinned_vertices; +} + +/* If the border is fixed, we simply pin the border vertices additionally to other pinned vertices. + */ +void GeometryData::retrieve_pinned_vertices(bool border_vertices_are_pinned) +{ + if (border_vertices_are_pinned) { + combine_matrices_of_pinned_and_boundary_vertices(); + } + else { + pinned_vertex_indices = VectorXi(explicitly_pinned_vertex_indices); + positions_of_pinned_vertices2d = MatrixXd(positions_of_explicitly_pinned_vertices2d); + } +} + +void GeometryData::initialize_if_needed(SLIMData &slim_data) const +{ + BLI_assert(slim_data.valid); + + if (!slim_data.skipInitialization) { + initialize_uvs(slim_data); + } +} + +void GeometryData::initialize_uvs(SLIMData &slim_data) const +{ + Eigen::MatrixXd uv_positions_of_boundary(boundary_vertex_indices.rows(), 2); + map_vertices_to_convex_border(uv_positions_of_boundary); + + bool all_vertices_on_boundary = (slim_data.V_o.rows() == uv_positions_of_boundary.rows()); + if (all_vertices_on_boundary) { + slim_data.V_o = uv_positions_of_boundary; + return; + } + + mvc(faces_by_vertexindices, + vertex_positions3d, + edges_by_vertexindices, + edge_lengths, + boundary_vertex_indices, + uv_positions_of_boundary, + slim_data.V_o); +} + +} // namespace slim diff --git a/intern/slim/intern/geometry_data_retrieval.h b/intern/slim/intern/geometry_data_retrieval.h new file mode 100644 index 00000000000..b56aa741fe3 --- /dev/null +++ b/intern/slim/intern/geometry_data_retrieval.h @@ -0,0 +1,70 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include + +#include + +#include "slim.h" +#include "slim_matrix_transfer.h" +#include "uv_initializer.h" + +using namespace Eigen; + +namespace slim { + +struct GeometryData { + int columns_2 = 2; + int columns_3 = 3; + int number_of_vertices = 0; + int number_of_faces = 0; + int number_of_edges_twice = 0; + int number_of_boundary_vertices = 0; + int number_of_pinned_vertices = 0; + + bool use_weights = false; + double weight_influence = 0.0; + + /* All the following maps have to be declared as last members. */ + Map vertex_positions3d = Map(NULL, 0, 0); + Map uv_positions2d = Map(NULL, 0, 0); + MatrixXd positions_of_pinned_vertices2d; + Map> positions_of_explicitly_pinned_vertices2d = + Map>(NULL, 0, 0); + + Map faces_by_vertexindices = Map(NULL, 0, 0); + Map edges_by_vertexindices = Map(NULL, 0, 0); + VectorXi pinned_vertex_indices; + Map explicitly_pinned_vertex_indices = Map(NULL, 0); + + Map edge_lengths = Map(NULL, 0); + Map boundary_vertex_indices = Map(NULL, 0); + Map weights_per_vertex = Map(NULL, 0); + + GeometryData(const MatrixTransfer &mt, MatrixTransferChart &chart); + GeometryData(const GeometryData &) = delete; + GeometryData &operator=(const GeometryData &) = delete; + + void construct_slim_data(SLIMData &slim_data, + bool skip_initialization, + int reflection_mode) const; + + void retrieve_pinned_vertices(bool border_vertices_are_pinned); + + private: + void set_geometry_data_matrices(SLIMData &slim_data) const; + bool has_valid_preinitialized_map() const; + bool can_initialization_be_skipped(bool skip_initialization) const; + void combine_matrices_of_pinned_and_boundary_vertices(); + void initialize_if_needed(SLIMData &slim_data) const; + void initialize_uvs(SLIMData &slim_data) const; +}; + +} // namespace slim diff --git a/intern/slim/intern/least_squares_relocator.cpp b/intern/slim/intern/least_squares_relocator.cpp new file mode 100644 index 00000000000..e56dc52046a --- /dev/null +++ b/intern/slim/intern/least_squares_relocator.cpp @@ -0,0 +1,240 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#include "least_squares_relocator.h" + +#include "slim.h" + +#include + +#include "BLI_assert.h" + +namespace slim { + +using namespace Eigen; + +static void apply_transformation(SLIMData &slim_data, Matrix2d &transformation_matrix) +{ + BLI_assert(slim_data.valid); + + for (int i = 0; i < slim_data.V_o.rows(); i++) { + slim_data.V_o.row(i) = transformation_matrix * slim_data.V_o.row(i).transpose(); + } +} + +static void apply_translation(SLIMData &slim_data, Vector2d &translation_vector) +{ + BLI_assert(slim_data.valid); + + for (int i = 0; i < slim_data.V_o.rows(); i++) { + slim_data.V_o.row(i) = translation_vector.transpose() + slim_data.V_o.row(i); + } +} + +static void retrieve_positions_of_pinned_vertices_in_initialization( + const MatrixXd &all_uv_positions_in_initialization, + const VectorXi &indices_of_pinned_vertices, + MatrixXd &position_of_pinned_vertices_in_initialization) +{ + int i = 0; + for (VectorXi::InnerIterator it(indices_of_pinned_vertices, 0); it; ++it, i++) { + int vertex_index = it.value(); + position_of_pinned_vertices_in_initialization.row(i) = all_uv_positions_in_initialization.row( + vertex_index); + } +} + +static void flip_input_geometry(SLIMData &slim_data) +{ + BLI_assert(slim_data.valid); + + VectorXi temp = slim_data.F.col(0); + slim_data.F.col(0) = slim_data.F.col(2); + slim_data.F.col(2) = temp; +} + +static void compute_centroid(const MatrixXd &point_cloud, Vector2d ¢roid) +{ + centroid << point_cloud.col(0).sum(), point_cloud.col(1).sum(); + centroid /= point_cloud.rows(); +} + +/* Finds scaling matrix: + * + * T = |a 0| + * |0 a| + * + * s.t. if to each point p in the inizialized map the following is applied + * + * T*p + * + * We get the closest scaling of the positions of the vertices in the initialized map to the pinned + * vertices in a least squares sense. We find them by solving + * + * argmin_{t} At = p + * + * i.e.: + * + * | x_1 | |u_1| + * | . | | . | + * | . | | . | + * | x_n | |u_n| + * | y_1 | * | a | = |v_1| + * | . | | . | + * | . | | . | + * | y_n | |v_n| + * + * `t` is of dimension `1 x 1` and `p` of dimension `2*numberOfPinnedVertices x 1` + * is the vector holding the uv positions of the pinned vertices. */ +static void compute_least_squares_scaling(MatrixXd centered_pins, + MatrixXd centered_initialized_pins, + Matrix2d &transformation_matrix) +{ + int number_of_pinned_vertices = centered_pins.rows(); + + MatrixXd a = MatrixXd::Zero(number_of_pinned_vertices * 2, 1); + a << centered_initialized_pins.col(0), centered_initialized_pins.col(1); + + VectorXd p(2 * number_of_pinned_vertices); + p << centered_pins.col(0), centered_pins.col(1); + + VectorXd t = a.colPivHouseholderQr().solve(p); + t(0) = abs(t(0)); + transformation_matrix << t(0), 0, 0, t(0); +} + +static void comput_least_squares_rotation_scale_only(SLIMData &slim_data, + Vector2d &translation_vector, + Matrix2d &transformation_matrix, + bool is_flip_allowed) +{ + BLI_assert(slim_data.valid); + + MatrixXd position_of_initialized_pins(slim_data.b.rows(), 2); + retrieve_positions_of_pinned_vertices_in_initialization( + slim_data.V_o, slim_data.b, position_of_initialized_pins); + + Vector2d centroid_of_initialized; + compute_centroid(position_of_initialized_pins, centroid_of_initialized); + + Vector2d centroid_of_pins; + compute_centroid(slim_data.bc, centroid_of_pins); + + MatrixXd centered_initialized_pins = position_of_initialized_pins.rowwise().operator-( + centroid_of_initialized.transpose()); + MatrixXd centeredpins = slim_data.bc.rowwise().operator-(centroid_of_pins.transpose()); + + MatrixXd s = centered_initialized_pins.transpose() * centeredpins; + + JacobiSVD svd(s, ComputeFullU | ComputeFullV); + + Matrix2d vu_t = svd.matrixV() * svd.matrixU().transpose(); + + Matrix2d singular_values = Matrix2d::Identity(); + + bool contains_reflection = vu_t.determinant() < 0; + if (contains_reflection) { + if (!is_flip_allowed) { + singular_values(1, 1) = vu_t.determinant(); + } + else { + flip_input_geometry(slim_data); + } + } + + compute_least_squares_scaling(centeredpins, centered_initialized_pins, transformation_matrix); + + transformation_matrix = transformation_matrix * svd.matrixV() * singular_values * + svd.matrixU().transpose(); + + translation_vector = centroid_of_pins - transformation_matrix * centroid_of_initialized; +} + +static void compute_transformation_matrix2_pins(const SLIMData &slim_data, + Matrix2d &transformation_matrix) +{ + BLI_assert(slim_data.valid); + + Vector2d pinned_position_difference_vector = slim_data.bc.row(0) - slim_data.bc.row(1); + Vector2d initialized_position_difference_vector = slim_data.V_o.row(slim_data.b(0)) - + slim_data.V_o.row(slim_data.b(1)); + + double scale = pinned_position_difference_vector.norm() / + initialized_position_difference_vector.norm(); + + pinned_position_difference_vector.normalize(); + initialized_position_difference_vector.normalize(); + + /* TODO: sometimes rotates in wrong direction. */ + double cos_angle = pinned_position_difference_vector.dot(initialized_position_difference_vector); + double sin_angle = sqrt(1 - pow(cos_angle, 2)); + + transformation_matrix << cos_angle, -sin_angle, sin_angle, cos_angle; + transformation_matrix = (Matrix2d::Identity() * scale) * transformation_matrix; +} + +static void compute_translation1_pin(const SLIMData &slim_data, Vector2d &translation_vector) +{ + BLI_assert(slim_data.valid); + translation_vector = slim_data.bc.row(0) - slim_data.V_o.row(slim_data.b(0)); +} + +static void transform_initialized_map(SLIMData &slim_data) +{ + BLI_assert(slim_data.valid); + Matrix2d transformation_matrix; + Vector2d translation_vector; + + int number_of_pinned_vertices = slim_data.b.rows(); + + switch (number_of_pinned_vertices) { + case 0: + return; + case 1: /* Only translation is needed with one pin. */ + compute_translation1_pin(slim_data, translation_vector); + apply_translation(slim_data, translation_vector); + break; + case 2: + compute_transformation_matrix2_pins(slim_data, transformation_matrix); + apply_transformation(slim_data, transformation_matrix); + compute_translation1_pin(slim_data, translation_vector); + apply_translation(slim_data, translation_vector); + break; + default: + + bool flip_allowed = slim_data.reflection_mode == 0; + + comput_least_squares_rotation_scale_only( + slim_data, translation_vector, transformation_matrix, flip_allowed); + + apply_transformation(slim_data, transformation_matrix); + apply_translation(slim_data, translation_vector); + + break; + } +} + +static bool is_translation_needed(const SLIMData &slim_data) +{ + BLI_assert(slim_data.valid); + bool pinned_vertices_exist = (slim_data.b.rows() > 0); + bool was_initialized = !slim_data.skipInitialization; + return was_initialized && pinned_vertices_exist; +} + +void transform_initialization_if_necessary(SLIMData &slim_data) +{ + BLI_assert(slim_data.valid); + + if (!is_translation_needed(slim_data)) { + return; + } + + transform_initialized_map(slim_data); +} +} // namespace slim diff --git a/intern/slim/intern/least_squares_relocator.h b/intern/slim/intern/least_squares_relocator.h new file mode 100644 index 00000000000..9d497925cc4 --- /dev/null +++ b/intern/slim/intern/least_squares_relocator.h @@ -0,0 +1,19 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include "slim.h" + +#include + +namespace slim { + +void transform_initialization_if_necessary(SLIMData &slimData); + +} diff --git a/intern/slim/intern/slim.cpp b/intern/slim/intern/slim.cpp new file mode 100644 index 00000000000..cf22cd420bb --- /dev/null +++ b/intern/slim/intern/slim.cpp @@ -0,0 +1,807 @@ +/* SPDX-FileCopyrightText: 2016 Michael Rabinovich + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#include "slim.h" +#include "doublearea.h" +#include "flip_avoiding_line_search.h" + +#include "BLI_assert.h" +#include "BLI_math_base.h" /* M_PI */ + +#include + +#include +#include +#include +#include + +namespace slim { + +/* GRAD + * G = grad(V,F) + * + * Compute the numerical gradient operator + * + * Inputs: + * V #vertices by 3 list of mesh vertex positions + * F #faces by 3 list of mesh face indices + * uniform #boolean (default false) - Use a uniform mesh instead of the vertices V + * Outputs: + * G #faces*dim by #V Gradient operator + * + * + * Gradient of a scalar function defined on piecewise linear elements (mesh) + * is constant on each triangle i,j,k: + * grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A + * where Xi is the scalar value at vertex i, Vi is the 3D position of vertex + * i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of + * 90 degrees + */ +template +static inline void grad(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::SparseMatrix &G, + bool uniform = false) +{ + Eigen::Matrix eperp21(F.rows(), 3), + eperp13(F.rows(), 3); + + for (int i = 0; i < F.rows(); ++i) { + /* Renaming indices of vertices of triangles for convenience. */ + int i1 = F(i, 0); + int i2 = F(i, 1); + int i3 = F(i, 2); + + /* #F x 3 matrices of triangle edge vectors, named after opposite vertices. */ + Eigen::Matrix v32 = V.row(i3) - V.row(i2); + Eigen::Matrix v13 = V.row(i1) - V.row(i3); + Eigen::Matrix v21 = V.row(i2) - V.row(i1); + Eigen::Matrix n = v32.cross(v13); + /* Area of parallelogram is twice area of triangle. + * Area of parallelogram is || v1 x v2 ||. + * This does correct l2 norm of rows, so that it contains #F list of twice. + * triangle areas. */ + double dblA = std::sqrt(n.dot(n)); + Eigen::Matrix u; + if (!uniform) { + /* Now normalize normals to get unit normals. */ + u = n / dblA; + } + else { + /* Abstract equilateral triangle v1=(0,0), v2=(h,0), v3=(h/2, (sqrt(3)/2)*h) */ + + /* Get h (by the area of the triangle). */ + double h = sqrt((dblA) / + sin(M_PI / 3.0)); /* (h^2*sin(60))/2. = Area => h = sqrt(2*Area/sin_60) */ + + Eigen::VectorXd v1, v2, v3; + v1 << 0, 0, 0; + v2 << h, 0, 0; + v3 << h / 2., (sqrt(3) / 2.) * h, 0; + + /* Now fix v32,v13,v21 and the normal. */ + v32 = v3 - v2; + v13 = v1 - v3; + v21 = v2 - v1; + n = v32.cross(v13); + } + + /* Rotate each vector 90 degrees around normal. */ + double norm21 = std::sqrt(v21.dot(v21)); + double norm13 = std::sqrt(v13.dot(v13)); + eperp21.row(i) = u.cross(v21); + eperp21.row(i) = eperp21.row(i) / std::sqrt(eperp21.row(i).dot(eperp21.row(i))); + eperp21.row(i) *= norm21 / dblA; + eperp13.row(i) = u.cross(v13); + eperp13.row(i) = eperp13.row(i) / std::sqrt(eperp13.row(i).dot(eperp13.row(i))); + eperp13.row(i) *= norm13 / dblA; + } + + std::vector rs; + rs.reserve(F.rows() * 4 * 3); + std::vector cs; + cs.reserve(F.rows() * 4 * 3); + std::vector vs; + vs.reserve(F.rows() * 4 * 3); + + /* Row indices. */ + for (int r = 0; r < 3; r++) { + for (int j = 0; j < 4; j++) { + for (int i = r * F.rows(); i < (r + 1) * F.rows(); i++) + rs.push_back(i); + } + } + + /* Column indices. */ + for (int r = 0; r < 3; r++) { + for (int i = 0; i < F.rows(); i++) + cs.push_back(F(i, 1)); + for (int i = 0; i < F.rows(); i++) + cs.push_back(F(i, 0)); + for (int i = 0; i < F.rows(); i++) + cs.push_back(F(i, 2)); + for (int i = 0; i < F.rows(); i++) + cs.push_back(F(i, 0)); + } + + /* Values. */ + for (int i = 0; i < F.rows(); i++) + vs.push_back(eperp13(i, 0)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(-eperp13(i, 0)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(eperp21(i, 0)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(-eperp21(i, 0)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(eperp13(i, 1)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(-eperp13(i, 1)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(eperp21(i, 1)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(-eperp21(i, 1)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(eperp13(i, 2)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(-eperp13(i, 2)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(eperp21(i, 2)); + for (int i = 0; i < F.rows(); i++) + vs.push_back(-eperp21(i, 2)); + + /* Create sparse gradient operator matrix.. */ + G.resize(3 * F.rows(), V.rows()); + std::vector> triplets; + for (int i = 0; i < (int)vs.size(); ++i) { + triplets.push_back(Eigen::Triplet(rs[i], cs[i], vs[i])); + } + G.setFromTriplets(triplets.begin(), triplets.end()); +} + +/* Computes the polar decomposition (R,T) of a matrix A using SVD singular + * value decomposition + * + * Inputs: + * A 3 by 3 matrix to be decomposed + * Outputs: + * R 3 by 3 rotation matrix part of decomposition (**always rotataion**) + * T 3 by 3 stretch matrix part of decomposition + * U 3 by 3 left-singular vectors + * S 3 by 1 singular values + * V 3 by 3 right-singular vectors + */ +template +static inline void polar_svd(const Eigen::PlainObjectBase &A, + Eigen::PlainObjectBase &R, + Eigen::PlainObjectBase &T, + Eigen::PlainObjectBase &U, + Eigen::PlainObjectBase &S, + Eigen::PlainObjectBase &V) +{ + using namespace std; + Eigen::JacobiSVD svd; + svd.compute(A, Eigen::ComputeFullU | Eigen::ComputeFullV); + U = svd.matrixU(); + V = svd.matrixV(); + S = svd.singularValues(); + R = U * V.transpose(); + const auto &SVT = S.asDiagonal() * V.adjoint(); + /* Check for reflection. */ + if (R.determinant() < 0) { + /* Annoyingly the .eval() is necessary. */ + auto W = V.eval(); + W.col(V.cols() - 1) *= -1.; + R = U * W.transpose(); + T = W * SVT; + } + else { + T = V * SVT; + } +} + +static inline void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, + const Eigen::MatrixXi &F, + const Eigen::MatrixXd &F1, + const Eigen::MatrixXd &F2, + Eigen::SparseMatrix &D1, + Eigen::SparseMatrix &D2) +{ + Eigen::SparseMatrix G; + grad(V, F, G); + Eigen::SparseMatrix Dx = G.block(0, 0, F.rows(), V.rows()); + Eigen::SparseMatrix Dy = G.block(F.rows(), 0, F.rows(), V.rows()); + Eigen::SparseMatrix Dz = G.block(2 * F.rows(), 0, F.rows(), V.rows()); + + D1 = F1.col(0).asDiagonal() * Dx + F1.col(1).asDiagonal() * Dy + F1.col(2).asDiagonal() * Dz; + D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy + F2.col(2).asDiagonal() * Dz; +} + +static inline void compute_weighted_jacobians(SLIMData &s, const Eigen::MatrixXd &uv) +{ + BLI_assert(s.valid); + + /* Ji=[D1*u,D2*u,D1*v,D2*v] */ + s.Ji.col(0) = s.Dx * uv.col(0); + s.Ji.col(1) = s.Dy * uv.col(0); + s.Ji.col(2) = s.Dx * uv.col(1); + s.Ji.col(3) = s.Dy * uv.col(1); + + /* Add weights. */ + Eigen::VectorXd weights = s.weightPerFaceMap.cast(); + s.Ji.col(0) = weights.cwiseProduct(s.Ji.col(0)); + s.Ji.col(1) = weights.cwiseProduct(s.Ji.col(1)); + s.Ji.col(2) = weights.cwiseProduct(s.Ji.col(2)); + s.Ji.col(3) = weights.cwiseProduct(s.Ji.col(3)); +} + +static inline void compute_unweighted_jacobians(SLIMData &s, const Eigen::MatrixXd &uv) +{ + BLI_assert(s.valid); + + /* Ji=[D1*u,D2*u,D1*v,D2*v] */ + s.Ji.col(0) = s.Dx * uv.col(0); + s.Ji.col(1) = s.Dy * uv.col(0); + s.Ji.col(2) = s.Dx * uv.col(1); + s.Ji.col(3) = s.Dy * uv.col(1); +} + +static inline void compute_jacobians(SLIMData &s, const Eigen::MatrixXd &uv) +{ + BLI_assert(s.valid); + + if (s.withWeightedParameterization) { + compute_weighted_jacobians(s, uv); + } + else { + compute_unweighted_jacobians(s, uv); + } +} + +static inline void update_weights_and_closest_rotations(SLIMData &s, Eigen::MatrixXd &uv) +{ + BLI_assert(s.valid); + compute_jacobians(s, uv); + + const double eps = 1e-8; + double exp_f = s.exp_factor; + + for (int i = 0; i < s.Ji.rows(); ++i) { + typedef Eigen::Matrix Mat2; + typedef Eigen::Matrix Vec2; + Mat2 ji, ri, ti, ui, vi; + Vec2 sing; + Vec2 closest_sing_vec; + Mat2 mat_W; + Vec2 m_sing_new; + double s1, s2; + + ji(0, 0) = s.Ji(i, 0); + ji(0, 1) = s.Ji(i, 1); + ji(1, 0) = s.Ji(i, 2); + ji(1, 1) = s.Ji(i, 3); + + polar_svd(ji, ri, ti, ui, sing, vi); + + s1 = sing(0); + s2 = sing(1); + + /* Update Weights according to energy. */ + switch (s.slim_energy) { + case SLIMData::ARAP: { + m_sing_new << 1, 1; + break; + } + case SLIMData::SYMMETRIC_DIRICHLET: { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + case SLIMData::LOG_ARAP: { + double s1_g = 2 * (log(s1) / s1); + double s2_g = 2 * (log(s2) / s2); + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + case SLIMData::CONFORMAL: { + double s1_g = 1 / (2 * s2) - s2 / (2 * pow(s1, 2)); + double s2_g = 1 / (2 * s1) - s1 / (2 * pow(s2, 2)); + + double geo_avg = sqrt(s1 * s2); + double s1_min = geo_avg; + double s2_min = geo_avg; + + m_sing_new << sqrt(s1_g / (2 * (s1 - s1_min))), sqrt(s2_g / (2 * (s2 - s2_min))); + + /* Change local step. */ + closest_sing_vec << s1_min, s2_min; + ri = ui * closest_sing_vec.asDiagonal() * vi.transpose(); + break; + } + case SLIMData::EXP_CONFORMAL: { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + + double in_exp = exp_f * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)); + double exp_thing = exp(in_exp); + + s1_g *= exp_thing * exp_f; + s2_g *= exp_thing * exp_f; + + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + case SLIMData::EXP_SYMMETRIC_DIRICHLET: { + double s1_g = 2 * (s1 - pow(s1, -3)); + double s2_g = 2 * (s2 - pow(s2, -3)); + + double in_exp = exp_f * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)); + double exp_thing = exp(in_exp); + + s1_g *= exp_thing * exp_f; + s2_g *= exp_thing * exp_f; + + m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1))); + break; + } + } + + if (std::abs(s1 - 1) < eps) + m_sing_new(0) = 1; + if (std::abs(s2 - 1) < eps) + m_sing_new(1) = 1; + mat_W = ui * m_sing_new.asDiagonal() * ui.transpose(); + + s.W_11(i) = mat_W(0, 0); + s.W_12(i) = mat_W(0, 1); + s.W_21(i) = mat_W(1, 0); + s.W_22(i) = mat_W(1, 1); + + /* 2) Update local step (doesn't have to be a rotation, for instance in case of conformal + * energy). */ + s.Ri(i, 0) = ri(0, 0); + s.Ri(i, 1) = ri(1, 0); + s.Ri(i, 2) = ri(0, 1); + s.Ri(i, 3) = ri(1, 1); + } +} + +template +static inline void local_basis(const Eigen::PlainObjectBase &V, + const Eigen::PlainObjectBase &F, + Eigen::PlainObjectBase &B1, + Eigen::PlainObjectBase &B2, + Eigen::PlainObjectBase &B3) +{ + using namespace Eigen; + using namespace std; + B1.resize(F.rows(), 3); + B2.resize(F.rows(), 3); + B3.resize(F.rows(), 3); + + for (unsigned i = 0; i < F.rows(); ++i) { + Eigen::Matrix v1 = + (V.row(F(i, 1)) - V.row(F(i, 0))).normalized(); + Eigen::Matrix t = V.row(F(i, 2)) - V.row(F(i, 0)); + Eigen::Matrix v3 = v1.cross(t).normalized(); + Eigen::Matrix v2 = v1.cross(v3).normalized(); + + B1.row(i) = v1; + B2.row(i) = -v2; + B3.row(i) = v3; + } +} + +static inline void pre_calc(SLIMData &s) +{ + BLI_assert(s.valid); + if (!s.has_pre_calc) { + s.v_n = s.v_num; + s.f_n = s.f_num; + + s.dim = 2; + Eigen::MatrixXd F1, F2, F3; + local_basis(s.V, s.F, F1, F2, F3); + compute_surface_gradient_matrix(s.V, s.F, F1, F2, s.Dx, s.Dy); + + s.W_11.resize(s.f_n); + s.W_12.resize(s.f_n); + s.W_21.resize(s.f_n); + s.W_22.resize(s.f_n); + + s.Dx.makeCompressed(); + s.Dy.makeCompressed(); + s.Dz.makeCompressed(); + s.Ri.resize(s.f_n, s.dim * s.dim); + s.Ji.resize(s.f_n, s.dim * s.dim); + s.rhs.resize(s.dim * s.v_num); + + /* Flattened weight matrix. */ + s.WGL_M.resize(s.dim * s.dim * s.f_n); + for (int i = 0; i < s.dim * s.dim; i++) + for (int j = 0; j < s.f_n; j++) + s.WGL_M(i * s.f_n + j) = s.M(j); + + s.first_solve = true; + s.has_pre_calc = true; + } +} + +static inline void buildA(SLIMData &s, Eigen::SparseMatrix &A) +{ + BLI_assert(s.valid); + /* Formula (35) in paper. */ + std::vector> IJV; + + IJV.reserve(4 * (s.Dx.outerSize() + s.Dy.outerSize())); + + /* A = [W11*Dx, W12*Dx; + * W11*Dy, W12*Dy; + * W21*Dx, W22*Dx; + * W21*Dy, W22*Dy]; */ + for (int k = 0; k < s.Dx.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(s.Dx, k); it; ++it) { + int dx_r = it.row(); + int dx_c = it.col(); + double val = it.value(); + double weight = s.weightPerFaceMap(dx_r); + + IJV.push_back(Eigen::Triplet(dx_r, dx_c, weight * val * s.W_11(dx_r))); + IJV.push_back(Eigen::Triplet(dx_r, s.v_n + dx_c, weight * val * s.W_12(dx_r))); + + IJV.push_back(Eigen::Triplet(2 * s.f_n + dx_r, dx_c, weight * val * s.W_21(dx_r))); + IJV.push_back( + Eigen::Triplet(2 * s.f_n + dx_r, s.v_n + dx_c, weight * val * s.W_22(dx_r))); + } + } + + for (int k = 0; k < s.Dy.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(s.Dy, k); it; ++it) { + int dy_r = it.row(); + int dy_c = it.col(); + double val = it.value(); + double weight = s.weightPerFaceMap(dy_r); + + IJV.push_back(Eigen::Triplet(s.f_n + dy_r, dy_c, weight * val * s.W_11(dy_r))); + IJV.push_back( + Eigen::Triplet(s.f_n + dy_r, s.v_n + dy_c, weight * val * s.W_12(dy_r))); + + IJV.push_back(Eigen::Triplet(3 * s.f_n + dy_r, dy_c, weight * val * s.W_21(dy_r))); + IJV.push_back( + Eigen::Triplet(3 * s.f_n + dy_r, s.v_n + dy_c, weight * val * s.W_22(dy_r))); + } + } + + A.setFromTriplets(IJV.begin(), IJV.end()); +} + +static inline void buildRhs(SLIMData &s, const Eigen::SparseMatrix &At) +{ + BLI_assert(s.valid); + + Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n); + f_rhs.setZero(); + + /* b = [W11*R11 + W12*R21; (formula (36)) + * W11*R12 + W12*R22; + * W21*R11 + W22*R21; + * W21*R12 + W22*R22]; */ + for (int i = 0; i < s.f_n; i++) { + f_rhs(i + 0 * s.f_n) = s.W_11(i) * s.Ri(i, 0) + s.W_12(i) * s.Ri(i, 1); + f_rhs(i + 1 * s.f_n) = s.W_11(i) * s.Ri(i, 2) + s.W_12(i) * s.Ri(i, 3); + f_rhs(i + 2 * s.f_n) = s.W_21(i) * s.Ri(i, 0) + s.W_22(i) * s.Ri(i, 1); + f_rhs(i + 3 * s.f_n) = s.W_21(i) * s.Ri(i, 2) + s.W_22(i) * s.Ri(i, 3); + } + + Eigen::VectorXd uv_flat(s.dim * s.v_n); + for (int i = 0; i < s.dim; i++) + for (int j = 0; j < s.v_n; j++) + uv_flat(s.v_n * i + j) = s.V_o(j, i); + + s.rhs = (At * s.WGL_M.asDiagonal() * f_rhs + s.proximal_p * uv_flat); +} + +static inline void add_soft_constraints(SLIMData &s, Eigen::SparseMatrix &L) +{ + BLI_assert(s.valid); + int v_n = s.v_num; + for (int d = 0; d < s.dim; d++) { + for (int i = 0; i < s.b.rows(); i++) { + int v_idx = s.b(i); + s.rhs(d * v_n + v_idx) += s.soft_const_p * s.bc(i, d); /* Right hand side. */ + L.coeffRef(d * v_n + v_idx, d * v_n + v_idx) += s.soft_const_p; /* Diagonal of matrix. */ + } + } +} + +static inline void build_linear_system(SLIMData &s, Eigen::SparseMatrix &L) +{ + BLI_assert(s.valid); + /* Formula (35) in paper. */ + Eigen::SparseMatrix A(s.dim * s.dim * s.f_n, s.dim * s.v_n); + buildA(s, A); + + Eigen::SparseMatrix At = A.transpose(); + At.makeCompressed(); + + Eigen::SparseMatrix id_m(At.rows(), At.rows()); + id_m.setIdentity(); + + /* Add proximal penalty. */ + L = At * s.WGL_M.asDiagonal() * A + s.proximal_p * id_m; /* Add also a proximal term. */ + L.makeCompressed(); + + buildRhs(s, At); + Eigen::SparseMatrix OldL = L; + add_soft_constraints(s, L); + L.makeCompressed(); +} + +static inline double compute_energy_with_jacobians(SLIMData &s, + const Eigen::MatrixXd &Ji, + Eigen::VectorXd &areas, + Eigen::VectorXd &singularValues, + bool gatherSingularValues) +{ + BLI_assert(s.valid); + double energy = 0; + + Eigen::Matrix ji; + for (int i = 0; i < s.f_n; i++) { + ji(0, 0) = Ji(i, 0); + ji(0, 1) = Ji(i, 1); + ji(1, 0) = Ji(i, 2); + ji(1, 1) = Ji(i, 3); + + typedef Eigen::Matrix Mat2; + typedef Eigen::Matrix Vec2; + Mat2 ri, ti, ui, vi; + Vec2 sing; + polar_svd(ji, ri, ti, ui, sing, vi); + double s1 = sing(0); + double s2 = sing(1); + + switch (s.slim_energy) { + case SLIMData::ARAP: { + energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2)); + break; + } + case SLIMData::SYMMETRIC_DIRICHLET: { + energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)); + + if (gatherSingularValues) { + singularValues(i) = s1; + singularValues(i + s.F.rows()) = s2; + } + break; + } + case SLIMData::EXP_SYMMETRIC_DIRICHLET: { + energy += areas(i) * + exp(s.exp_factor * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2))); + break; + } + case SLIMData::LOG_ARAP: { + energy += areas(i) * (pow(log(s1), 2) + pow(log(s2), 2)); + break; + } + case SLIMData::CONFORMAL: { + energy += areas(i) * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)); + break; + } + case SLIMData::EXP_CONFORMAL: { + energy += areas(i) * exp(s.exp_factor * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2))); + break; + } + } + } + + return energy; +} + +static inline double compute_soft_const_energy(SLIMData &s, Eigen::MatrixXd &V_o) +{ + BLI_assert(s.valid); + double e = 0; + for (int i = 0; i < s.b.rows(); i++) { + e += s.soft_const_p * (s.bc.row(i) - V_o.row(s.b(i))).squaredNorm(); + } + return e; +} + +static inline double compute_energy(SLIMData &s, + Eigen::MatrixXd &V_new, + Eigen::VectorXd &singularValues, + bool gatherSingularValues) +{ + BLI_assert(s.valid); + compute_jacobians(s, V_new); + return compute_energy_with_jacobians(s, s.Ji, s.M, singularValues, gatherSingularValues) + + compute_soft_const_energy(s, V_new); +} + +static inline double compute_energy(SLIMData &s, Eigen::MatrixXd &V_new) +{ + BLI_assert(s.valid); + Eigen::VectorXd temp; + return compute_energy(s, V_new, temp, false); +} + +static inline double compute_energy(SLIMData &s, + Eigen::MatrixXd &V_new, + Eigen::VectorXd &singularValues) +{ + BLI_assert(s.valid); + return compute_energy(s, V_new, singularValues, true); +} + +void slim_precompute(Eigen::MatrixXd &V, + Eigen::MatrixXi &F, + Eigen::MatrixXd &V_init, + SLIMData &data, + SLIMData::SLIM_ENERGY slim_energy, + Eigen::VectorXi &b, + Eigen::MatrixXd &bc, + double soft_p) +{ + BLI_assert(data.valid); + data.V = V; + data.F = F; + data.V_o = V_init; + + data.v_num = V.rows(); + data.f_num = F.rows(); + + data.slim_energy = slim_energy; + + data.b = b; + data.bc = bc; + data.soft_const_p = soft_p; + + data.proximal_p = 0.0001; + + doublearea(V, F, data.M); + data.M /= 2.; + data.mesh_area = data.M.sum(); + + data.mesh_improvement_3d = false; /* Whether to use a jacobian derived from a real mesh or an + * abstract regular mesh (used for mesh improvement). */ + data.exp_factor = + 1.0; /* Param used only for exponential energies (e.g exponential symmetric dirichlet). */ + + assert(F.cols() == 3); + + pre_calc(data); + + data.energy = compute_energy(data, data.V_o) / data.mesh_area; +} + +inline double computeGlobalScaleInvarianceFactor(Eigen::VectorXd &singularValues, + Eigen::VectorXd &areas) +{ + int nFaces = singularValues.rows() / 2; + + Eigen::VectorXd areasChained(2 * nFaces); + areasChained << areas, areas; + + /* Per face energy for face i with singvals si1 and si2 and area ai when scaling geometry by x is + * + * ai*(si1*x)^2 + ai*(si2*x)^2 + ai/(si1*x)^2 + ai/(si2*x)^2) + * + * The combined Energy of all faces is therefore + * (s1 and s2 are the sums over all ai*(si1^2) and ai*(si2^2) respectively. t1 and t2 + * are the sums over all ai/(si1^2) and ai/(si2^2) respectively) + * + * s1*(x^2) + s2*(x^2) + t1/(x^2) + t2/(x^2) + * + * with a = (s1 + s2) and b = (t1 + t2) we get + * + * ax^2 + b/x^2 + * + * it's derivative is + * + * 2ax - 2b/(x^3) + * + * and when we set it zero we get + * + * x^4 = b/a => x = sqrt(sqrt(b/a)) + */ + + Eigen::VectorXd squaredSingularValues = singularValues.cwiseProduct(singularValues); + Eigen::VectorXd inverseSquaredSingularValues = + singularValues.cwiseProduct(singularValues).cwiseInverse(); + + Eigen::VectorXd weightedSquaredSingularValues = squaredSingularValues.cwiseProduct(areasChained); + Eigen::VectorXd weightedInverseSquaredSingularValues = inverseSquaredSingularValues.cwiseProduct( + areasChained); + + double s1 = weightedSquaredSingularValues.head(nFaces).sum(); + double s2 = weightedSquaredSingularValues.tail(nFaces).sum(); + + double t1 = weightedInverseSquaredSingularValues.head(nFaces).sum(); + double t2 = weightedInverseSquaredSingularValues.tail(nFaces).sum(); + + double a = s1 + s2; + double b = t1 + t2; + + double x = sqrt(sqrt(b / a)); + + return 1 / x; +} + +static inline void solve_weighted_arap(SLIMData &s, Eigen::MatrixXd &uv) +{ + BLI_assert(s.valid); + using namespace Eigen; + + Eigen::SparseMatrix L; + build_linear_system(s, L); + + /* Solve. */ + Eigen::VectorXd Uc; + SimplicialLDLT> solver; + Uc = solver.compute(L).solve(s.rhs); + + for (int i = 0; i < Uc.size(); i++) + if (!std::isfinite(Uc(i))) + throw SlimFailedException(); + + for (int i = 0; i < s.dim; i++) + uv.col(i) = Uc.block(i * s.v_n, 0, s.v_n, 1); +} + +Eigen::MatrixXd slim_solve(SLIMData &data, int iter_num) +{ + BLI_assert(data.valid); + Eigen::VectorXd singularValues; + bool are_pins_present = data.b.rows() > 0; + + if (are_pins_present) { + singularValues.resize(data.F.rows() * 2); + data.energy = compute_energy(data, data.V_o, singularValues) / data.mesh_area; + } + + for (int i = 0; i < iter_num; i++) { + Eigen::MatrixXd dest_res; + dest_res = data.V_o; + + /* Solve Weighted Proxy. */ + update_weights_and_closest_rotations(data, dest_res); + solve_weighted_arap(data, dest_res); + + std::function compute_energy_func = [&](Eigen::MatrixXd &aaa) { + return are_pins_present ? compute_energy(data, aaa, singularValues) : + compute_energy(data, aaa); + }; + + data.energy = flip_avoiding_line_search(data.F, + data.V_o, + dest_res, + compute_energy_func, + data.energy * data.mesh_area) / + data.mesh_area; + + if (are_pins_present) { + data.globalScaleInvarianceFactor = computeGlobalScaleInvarianceFactor(singularValues, + data.M); + data.Dx /= data.globalScaleInvarianceFactor; + data.Dy /= data.globalScaleInvarianceFactor; + data.energy = compute_energy(data, data.V_o, singularValues) / data.mesh_area; + } + } + + return data.V_o; +} + +} // namespace slim diff --git a/intern/slim/intern/slim.h b/intern/slim/intern/slim.h new file mode 100644 index 00000000000..1897745008f --- /dev/null +++ b/intern/slim/intern/slim.h @@ -0,0 +1,122 @@ +/* SPDX-FileCopyrightText: 2016 Michael Rabinovich + * 2023 Blender Authors + * + * SPDX-License-Identifier: MPL-2.0 */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include +#include + +#include + +namespace slim { + +class SlimFailedException : public std::runtime_error { + public: + SlimFailedException() : std::runtime_error("Slim operation failed") {} +}; + +/* Compute a SLIM map as derived in "Scalable Locally Injective Maps" [Rabinovich et al. 2016].. */ +struct SLIMData { + bool valid = true; + + /* Input. */ + Eigen::MatrixXd V; /* #V by 3 list of mesh vertex positions. */ + Eigen::MatrixXi F; /* #F by 3/3 list of mesh faces (triangles). */ + enum SLIM_ENERGY { + ARAP, + LOG_ARAP, + SYMMETRIC_DIRICHLET, + CONFORMAL, + EXP_CONFORMAL, + EXP_SYMMETRIC_DIRICHLET + }; + SLIM_ENERGY slim_energy; + + /* Optional Input. */ + /* Soft constraints. */ + Eigen::VectorXi b; + Eigen::MatrixXd bc; + double soft_const_p; + + double exp_factor; /* Used for exponential energies, ignored otherwise. */ + bool mesh_improvement_3d; /* Only supported for 3d. */ + + int reflection_mode; + bool skipInitialization = false; + bool validPreInitialization = false; + double expectedSurfaceAreaOfResultingMap = 0; + + /* Output. */ + Eigen::MatrixXd V_o; /* #V by dim list of mesh vertex positions (dim = 2 for parametrization, 3 + otherwise). */ + Eigen::MatrixXd + oldUVs; /* #V by dim list of mesh vertex positions (dim = 2 for parametrization,. */ + /* 3 otherwise). */ + + /* Weight-map for weighted parameterization. */ + bool withWeightedParameterization; + Eigen::VectorXf weightmap; + Eigen::VectorXf weightPerFaceMap; + double weightInfluence; + double globalScaleInvarianceFactor = 1.0; + + double energy; /* Objective value. */ + + /* Internal. */ + Eigen::VectorXd M; + double mesh_area; + double avg_edge_length; + int v_num; + int f_num; + double proximal_p; + + Eigen::VectorXd WGL_M; + Eigen::VectorXd rhs; + Eigen::MatrixXd Ri, Ji; + Eigen::VectorXd W_11; + Eigen::VectorXd W_12; + Eigen::VectorXd W_13; + Eigen::VectorXd W_21; + Eigen::VectorXd W_22; + Eigen::VectorXd W_23; + Eigen::VectorXd W_31; + Eigen::VectorXd W_32; + Eigen::VectorXd W_33; + Eigen::SparseMatrix Dx, Dy, Dz; + int f_n, v_n; + bool first_solve; + bool has_pre_calc = false; + int dim; +}; + +/* Compute necessary information to start using SLIM + * Inputs: + * V #V by 3 list of mesh vertex positions + * F #F by 3/3 list of mesh faces (triangles) + * b list of boundary indices into V + * bc #b by dim list of boundary conditions + * soft_p Soft penalty factor (can be zero) + * slim_energy Energy to minimize + */ +void slim_precompute(Eigen::MatrixXd &V, + Eigen::MatrixXi &F, + Eigen::MatrixXd &V_init, + SLIMData &data, + SLIMData::SLIM_ENERGY slim_energy, + Eigen::VectorXi &b, + Eigen::MatrixXd &bc, + double soft_p); + +/* Run iter_num iterations of SLIM + * Outputs: + * V_o (in SLIMData): #V by dim list of mesh vertex positions + */ +Eigen::MatrixXd slim_solve(SLIMData &data, int iter_num); + +} // namespace slim diff --git a/intern/slim/intern/slim_matrix_transfer.cpp b/intern/slim/intern/slim_matrix_transfer.cpp new file mode 100644 index 00000000000..276e7c383db --- /dev/null +++ b/intern/slim/intern/slim_matrix_transfer.cpp @@ -0,0 +1,50 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#include + +#include "geometry_data_retrieval.h" +#include "slim.h" +#include "slim_matrix_transfer.h" + +namespace slim { + +MatrixTransferChart::MatrixTransferChart() = default; +MatrixTransferChart::MatrixTransferChart(MatrixTransferChart &&) = default; +MatrixTransferChart::~MatrixTransferChart() = default; +MatrixTransfer::MatrixTransfer() = default; +MatrixTransfer::~MatrixTransfer() = default; + +void MatrixTransferChart::free_slim_data() +{ + data.reset(nullptr); +} + +void MatrixTransfer::setup_slim_data(MatrixTransferChart &chart) const +{ + SLIMDataPtr slim_data = std::make_unique(); + + try { + if (!chart.succeeded) { + throw SlimFailedException(); + } + + GeometryData geometry_data(*this, chart); + geometry_data.construct_slim_data(*slim_data, skip_initialization, reflection_mode); + + chart.pinned_vertices_num = geometry_data.number_of_pinned_vertices; + } + catch (SlimFailedException &) { + slim_data->valid = false; + chart.succeeded = false; + } + + chart.data = std::move(slim_data); +} + +} // namespace slim diff --git a/intern/slim/intern/slim_parametrizer.cpp b/intern/slim/intern/slim_parametrizer.cpp new file mode 100644 index 00000000000..ea54337120f --- /dev/null +++ b/intern/slim/intern/slim_parametrizer.cpp @@ -0,0 +1,161 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#include +#include +#include + +#include "area_compensation.h" +#include "doublearea.h" +#include "geometry_data_retrieval.h" +#include "least_squares_relocator.h" +#include "slim.h" +#include "uv_initializer.h" + +#include +#include +#include + +namespace slim { + +using namespace Eigen; + +static void transfer_uvs_back_to_native_part(MatrixTransferChart &chart, Eigen::MatrixXd &uv) +{ + if (!chart.succeeded) { + return; + } + + auto &uv_coordinate_array = chart.uv_matrices; + int number_of_vertices = chart.verts_num; + + for (int i = 0; i < number_of_vertices; i++) { + uv_coordinate_array[i] = uv(i, 0); + uv_coordinate_array[number_of_vertices + i] = uv(i, 1); + } +} + +static Eigen::MatrixXd get_interactive_result_blended_with_original(float blend, + const SLIMData &slim_data) +{ + Eigen::MatrixXd original_map_weighted = blend * slim_data.oldUVs; + Eigen::MatrixXd interactive_result_map = (1.0 - blend) * slim_data.V_o; + return original_map_weighted + interactive_result_map; +} + +static void adjust_pins(SLIMData &slim_data, const PinnedVertexData &pinned_vertex_data) +{ + if (!slim_data.valid) { + return; + } + + auto &pinned_vertex_indices = pinned_vertex_data.pinned_vertex_indices; + auto &pinned_vertex_positions_2D = pinned_vertex_data.pinned_vertex_positions_2D; + auto &selected_pins = pinned_vertex_data.selected_pins; + + int n_pins = pinned_vertex_indices.size(); + int n_selected_pins = selected_pins.size(); + + Eigen::VectorXi old_pin_indices = slim_data.b; + Eigen::MatrixXd old_pin_positions = slim_data.bc; + + slim_data.b.resize(n_pins); + slim_data.bc.resize(n_pins, 2); + + int old_pin_pointer = 0; + int selected_pin_pointer = 0; + + for (int new_pin_pointer = 0; new_pin_pointer < n_pins; new_pin_pointer++) { + + int pinned_vertex_index = pinned_vertex_indices[new_pin_pointer]; + slim_data.b(new_pin_pointer) = pinned_vertex_index; + + while ((old_pin_pointer < old_pin_indices.size()) && + (old_pin_indices(old_pin_pointer) < pinned_vertex_index)) + { + ++old_pin_pointer; + } + bool old_pointer_valid = (old_pin_pointer < old_pin_indices.size()) && + (old_pin_indices(old_pin_pointer) == pinned_vertex_index); + + while ((selected_pin_pointer < n_selected_pins) && + (selected_pins[selected_pin_pointer] < pinned_vertex_index)) + { + ++selected_pin_pointer; + } + bool pin_selected = (selected_pin_pointer < n_selected_pins) && + (selected_pins[selected_pin_pointer] == pinned_vertex_index); + + if (!pin_selected && old_pointer_valid) { + slim_data.bc.row(new_pin_pointer) = old_pin_positions.row(old_pin_pointer); + } + else { + slim_data.bc(new_pin_pointer, 0) = pinned_vertex_positions_2D[2 * new_pin_pointer]; + slim_data.bc(new_pin_pointer, 1) = pinned_vertex_positions_2D[2 * new_pin_pointer + 1]; + } + } +} + +void MatrixTransferChart::transfer_uvs_blended(float blend) +{ + if (!succeeded) { + return; + } + + Eigen::MatrixXd blended_uvs = get_interactive_result_blended_with_original(blend, *data); + correct_map_surface_area_if_necessary(*data); + transfer_uvs_back_to_native_part(*this, blended_uvs); +} + +void MatrixTransferChart::try_slim_solve(int iter_num) +{ + if (!succeeded) { + return; + } + + try { + slim_solve(*data, iter_num); + } + catch (SlimFailedException &) { + succeeded = false; + } +} + +void MatrixTransferChart::parametrize_single_iteration() +{ + int number_of_iterations = 1; + try_slim_solve(number_of_iterations); +} + +void MatrixTransfer::parametrize_live(MatrixTransferChart &chart, + const PinnedVertexData &pinned_vertex_data) +{ + int number_of_iterations = 3; + adjust_pins(*chart.data, pinned_vertex_data); + + chart.try_slim_solve(number_of_iterations); + + correct_map_surface_area_if_necessary(*chart.data); + transfer_uvs_back_to_native_part(chart, chart.data->V_o); +} + +void MatrixTransfer::parametrize() +{ + for (MatrixTransferChart &chart : charts) { + setup_slim_data(chart); + + chart.try_slim_solve(n_iterations); + + correct_map_surface_area_if_necessary(*chart.data); + transfer_uvs_back_to_native_part(chart, chart.data->V_o); + + chart.free_slim_data(); + } +} + +} // namespace slim diff --git a/intern/slim/intern/uv_initializer.cpp b/intern/slim/intern/uv_initializer.cpp new file mode 100644 index 00000000000..233746b7bd9 --- /dev/null +++ b/intern/slim/intern/uv_initializer.cpp @@ -0,0 +1,293 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#include "uv_initializer.h" +#include "cotmatrix.h" + +#include + +namespace slim { + +static double compute_angle(const Eigen::Vector3d &a, const Eigen::Vector3d &b) +{ + return acos(a.dot(b) / (a.norm() * b.norm())); +} + +static void find_vertex_to_opposite_angles_correspondence( + const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + Eigen::SparseMatrix &vertex_to_face_indices) +{ + + typedef Eigen::Triplet t; + std::vector coefficients; + + for (int i = 0; i < f.rows(); i++) { + + int vertex_index1 = f(i, 0); + int vertex_index2 = f(i, 1); + int vertex_index3 = f(i, 2); + + double angle1 = compute_angle(v.row(vertex_index2) - v.row(vertex_index1), + v.row(vertex_index3) - v.row(vertex_index1)); + double angle2 = compute_angle(v.row(vertex_index3) - v.row(vertex_index2), + v.row(vertex_index1) - v.row(vertex_index2)); + double angle3 = compute_angle(v.row(vertex_index1) - v.row(vertex_index3), + v.row(vertex_index2) - v.row(vertex_index3)); + + coefficients.push_back(t(vertex_index1, 2 * vertex_index2, angle3)); + coefficients.push_back(t(vertex_index1, 2 * vertex_index3 + 1, angle2)); + + coefficients.push_back(t(vertex_index2, 2 * vertex_index1 + 1, angle3)); + coefficients.push_back(t(vertex_index2, 2 * vertex_index3, angle1)); + + coefficients.push_back(t(vertex_index3, 2 * vertex_index1, angle2)); + coefficients.push_back(t(vertex_index3, 2 * vertex_index2 + 1, angle1)); + } + + vertex_to_face_indices.setFromTriplets(coefficients.begin(), coefficients.end()); +} + +static void find_vertex_to_its_angles_correspondence( + const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + Eigen::SparseMatrix &vertex_to_face_indices) +{ + + typedef Eigen::Triplet t; + std::vector coefficients; + + for (int i = 0; i < f.rows(); i++) { + + int vertex_index1 = f(i, 0); + int vertex_index2 = f(i, 1); + int vertex_index3 = f(i, 2); + + double angle1 = compute_angle(v.row(vertex_index2) - v.row(vertex_index1), + v.row(vertex_index3) - v.row(vertex_index1)); + double angle2 = compute_angle(v.row(vertex_index3) - v.row(vertex_index2), + v.row(vertex_index1) - v.row(vertex_index2)); + double angle3 = compute_angle(v.row(vertex_index1) - v.row(vertex_index3), + v.row(vertex_index2) - v.row(vertex_index3)); + + coefficients.push_back(t(vertex_index1, 2 * vertex_index2, angle1)); + coefficients.push_back(t(vertex_index1, 2 * vertex_index3 + 1, angle1)); + + coefficients.push_back(t(vertex_index2, 2 * vertex_index1 + 1, angle2)); + coefficients.push_back(t(vertex_index2, 2 * vertex_index3, angle2)); + + coefficients.push_back(t(vertex_index3, 2 * vertex_index1, angle3)); + coefficients.push_back(t(vertex_index3, 2 * vertex_index2 + 1, angle3)); + } + + vertex_to_face_indices.setFromTriplets(coefficients.begin(), coefficients.end()); +} + +/* Implementation of different fixed-border parameterizations, mean value coordinates, harmonic, + * tutte. */ +void convex_border_parameterization(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &uv, + Method method) +{ + int verts_num = uv.rows(); + int edges_num = e.rows(); + + Eigen::SparseMatrix vertex_to_angles(verts_num, verts_num * 2); + + switch (method) { + case HARMONIC: + find_vertex_to_opposite_angles_correspondence(f, v, vertex_to_angles); + break; + case MVC: + find_vertex_to_its_angles_correspondence(f, v, vertex_to_angles); + break; + case TUTTE: + break; + } + + int n_unknowns = verts_num - bnd.size(); + int n_knowns = bnd.size(); + + Eigen::SparseMatrix aint(n_unknowns, n_unknowns); + Eigen::SparseMatrix abnd(n_unknowns, n_knowns); + Eigen::VectorXd z(n_knowns); + + std::vector> int_triplet_vector; + std::vector> bnd_triplet_vector; + + int rowindex; + int columnindex; + double edge_weight, edge_length; + Eigen::RowVector2i edge; + + int first_vertex, second_vertex; + + for (int e_idx = 0; e_idx < edges_num; e_idx++) { + edge = e.row(e_idx); + edge_length = el(e_idx); + first_vertex = edge(0); + second_vertex = edge(1); + + if (first_vertex >= n_knowns) { + /* Into aint. */ + rowindex = first_vertex - n_knowns; + + double angle1 = vertex_to_angles.coeff(first_vertex, 2 * second_vertex); + double angle2 = vertex_to_angles.coeff(first_vertex, 2 * second_vertex + 1); + + switch (method) { + case HARMONIC: + edge_weight = 1 / tan(angle1) + 1 / tan(angle2); + break; + case MVC: + edge_weight = tan(angle1 / 2) + tan(angle2 / 2); + edge_weight /= edge_length; + break; + case TUTTE: + edge_weight = 1; + break; + } + + int_triplet_vector.push_back(Eigen::Triplet(rowindex, rowindex, edge_weight)); + + if (second_vertex >= n_knowns) { + /* Also an unknown point in the interior. */ + columnindex = second_vertex - n_knowns; + + int_triplet_vector.push_back(Eigen::Triplet(rowindex, columnindex, -edge_weight)); + } + else { + /* Known point on the border. */ + columnindex = second_vertex; + bnd_triplet_vector.push_back(Eigen::Triplet(rowindex, columnindex, edge_weight)); + } + } + } + + aint.setFromTriplets(int_triplet_vector.begin(), int_triplet_vector.end()); + aint.makeCompressed(); + + abnd.setFromTriplets(bnd_triplet_vector.begin(), bnd_triplet_vector.end()); + abnd.makeCompressed(); + + for (int i = 0; i < n_unknowns; i++) { + double factor = aint.coeff(i, i); + aint.row(i) /= factor; + abnd.row(i) /= factor; + } + + Eigen::SparseLU> solver; + solver.compute(aint); + + for (int i = 0; i < 2; i++) { + + for (int zindex = 0; zindex < n_knowns; zindex++) { + z(zindex) = bnd_uv(bnd(zindex), i); + } + + Eigen::VectorXd b = abnd * z; + + Eigen::VectorXd uvs; + uvs = solver.solve(b); + + Eigen::VectorXd boundary = bnd_uv.col(i); + Eigen::VectorXd interior = uvs; + + uv.col(i) << boundary, interior; + } +} + +void mvc(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &uv) +{ + + convex_border_parameterization(f, v, e, el, bnd, bnd_uv, uv, Method::MVC); +} + +void harmonic(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &uv) +{ + + convex_border_parameterization(f, v, e, el, bnd, bnd_uv, uv, Method::HARMONIC); +} + +void tutte(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &uv) +{ + + convex_border_parameterization(f, v, e, el, bnd, bnd_uv, uv, Method::TUTTE); +} + +void map_vertices_to_convex_border(Eigen::MatrixXd &vertex_positions) +{ + double pi = atan(1) * 4; + int boundary_vertices_num = vertex_positions.rows(); + double x, y; + double angle = 2 * pi / boundary_vertices_num; + + for (int i = 0; i < boundary_vertices_num; i++) { + x = cos(angle * i); + y = sin(angle * i); + vertex_positions(i, 0) = (x * 0.5) + 0.5; + vertex_positions(i, 1) = (y * 0.5) + 0.5; + } +} + +static void get_flips(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &uv, + std::vector &flip_idx) +{ + flip_idx.resize(0); + for (int i = 0; i < f.rows(); i++) { + + Eigen::Vector2d v1_n = uv.row(f(i, 0)); + Eigen::Vector2d v2_n = uv.row(f(i, 1)); + Eigen::Vector2d v3_n = uv.row(f(i, 2)); + + Eigen::MatrixXd t2_homo(3, 3); + t2_homo.col(0) << v1_n(0), v1_n(1), 1; + t2_homo.col(1) << v2_n(0), v2_n(1), 1; + t2_homo.col(2) << v3_n(0), v3_n(1), 1; + double det = t2_homo.determinant(); + assert(det == det); + if (det < 0) { + flip_idx.push_back(i); + } + } +} + +int count_flips(const Eigen::MatrixXi &f, const Eigen::MatrixXd &uv) +{ + + std::vector flip_idx; + get_flips(f, uv, flip_idx); + + return flip_idx.size(); +} + +} // namespace slim diff --git a/intern/slim/intern/uv_initializer.h b/intern/slim/intern/uv_initializer.h new file mode 100644 index 00000000000..281943a7c09 --- /dev/null +++ b/intern/slim/intern/uv_initializer.h @@ -0,0 +1,64 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include + +#include +#include + +enum Method { TUTTE, HARMONIC, MVC }; + +namespace slim { + +void convex_border_parameterization(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &UV, + Method method); + +void mvc(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &UV); + +void harmonic(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &UV); + +void tutte(const Eigen::MatrixXi &f, + const Eigen::MatrixXd &v, + const Eigen::MatrixXi &e, + const Eigen::VectorXd &el, + const Eigen::VectorXi &bnd, + const Eigen::MatrixXd &bnd_uv, + Eigen::MatrixXd &UV); + +void harmonic(const Eigen::MatrixXd &v, + const Eigen::MatrixXi &f, + const Eigen::MatrixXi &B, + const Eigen::MatrixXd &bnd_uv, + int power_of_harmonic_operaton, + Eigen::MatrixXd &UV); + +void map_vertices_to_convex_border(Eigen::MatrixXd &vertex_positions); + +int count_flips(const Eigen::MatrixXi &f, const Eigen::MatrixXd &uv); + +} // namespace slim diff --git a/intern/slim/slim_matrix_transfer.h b/intern/slim/slim_matrix_transfer.h new file mode 100644 index 00000000000..2db517412e7 --- /dev/null +++ b/intern/slim/slim_matrix_transfer.h @@ -0,0 +1,110 @@ +/* SPDX-FileCopyrightText: 2023 Blender Authors + * + * SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup intern_slim + */ + +#pragma once + +#include +#include + +namespace slim { + +struct SLIMData; + +typedef std::unique_ptr SLIMDataPtr; + +/** + * MatrixTransferChart holds all information and data matrices to be + * transferred from Blender to SLIM. + */ +struct MatrixTransferChart { + int verts_num = 0; + int faces_num = 0; + int pinned_vertices_num = 0; + int boundary_vertices_num = 0; + int edges_num = 0; + + /** Field indicating whether a given SLIM operation succeeded or not. */ + bool succeeded = false; + + /** Vertex positions (matrix [verts_num x 3] in a linearized form). */ + std::vector v_matrices; + /** UV positions of vertices (matrix [verts_num x 2] in a linearized form). */ + std::vector uv_matrices; + /** Positions of pinned vertices (matrix [pinned_vertices_num x 2] in a linearized form). */ + std::vector pp_matrices; + /** Edge lengths. */ + std::vector el_vectors; + /** Weights per vertex. */ + std::vector w_vectors; + + /** Vertex index triplets making up faces (matrix [faces_num x 3] in a linearized form). */ + std::vector f_matrices; + /** Indices of pinned vertices. */ + std::vector p_matrices; + /** Vertex index tuples making up edges (matrix [edges_num x 2] in a linearized form). */ + std::vector e_matrices; + /** Vertex indices of boundary vertices. */ + std::vector b_vectors; + + SLIMDataPtr data; + + MatrixTransferChart(); + MatrixTransferChart(MatrixTransferChart &&); + + MatrixTransferChart(const MatrixTransferChart &) = delete; + MatrixTransferChart &operator=(const MatrixTransferChart &) = delete; + + ~MatrixTransferChart(); + + void try_slim_solve(int iter_num); + /** Executes a single iteration of SLIM, must follow a proper setup & initialization. */ + void parametrize_single_iteration(); + /** + * Called from the native part during each iteration of interactive parametrization. + * The blend parameter decides the linear blending between the original UV map and the one + * obtained from the accumulated SLIM iterations so far. + */ + void transfer_uvs_blended(float blend); + void free_slim_data(); +}; + +struct PinnedVertexData { + std::vector pinned_vertex_indices; + std::vector pinned_vertex_positions_2D; + std::vector selected_pins; +}; + +struct MatrixTransfer { + bool fixed_boundary = false; + bool use_weights = false; + double weight_influence = 0.0; + int reflection_mode = 0; + int n_iterations = 0; + bool skip_initialization = false; + bool is_minimize_stretch = false; + + std::vector charts; + + /** Used for pins update in live unwrap. */ + PinnedVertexData pinned_vertex_data; + + MatrixTransfer(); + MatrixTransfer(const MatrixTransfer &) = delete; + MatrixTransfer &operator=(const MatrixTransfer &) = delete; + ~MatrixTransfer(); + + void parametrize(); + + /** Executes slim iterations during live unwrap. needs to provide new selected-pin positions. */ + void parametrize_live(MatrixTransferChart &chart, const PinnedVertexData &pinned_vertex_data); + + /** Transfers all the matrices from the native part and initializes SLIM. */ + void setup_slim_data(MatrixTransferChart &chart) const; +}; + +} // namespace slim diff --git a/scripts/startup/bl_ui/space_image.py b/scripts/startup/bl_ui/space_image.py index 81007fc2289..f1c79169c23 100644 --- a/scripts/startup/bl_ui/space_image.py +++ b/scripts/startup/bl_ui/space_image.py @@ -402,7 +402,7 @@ class IMAGE_MT_uvs_unwrap(Menu): def draw(self, _context): layout = self.layout - layout.operator("uv.unwrap") + layout.operator_enum("uv.unwrap", "method") layout.separator() diff --git a/scripts/startup/bl_ui/space_view3d.py b/scripts/startup/bl_ui/space_view3d.py index 7b151798fa1..2b0a6081f44 100644 --- a/scripts/startup/bl_ui/space_view3d.py +++ b/scripts/startup/bl_ui/space_view3d.py @@ -1595,21 +1595,7 @@ class VIEW3D_MT_uv_map(Menu): def draw(self, _context): layout = self.layout - layout.operator("uv.unwrap") - - layout.separator() - - layout.operator_context = 'INVOKE_DEFAULT' - layout.operator("uv.smart_project") - layout.operator("uv.lightmap_pack") - layout.operator("uv.follow_active_quads") - - layout.separator() - - layout.operator_context = 'EXEC_REGION_WIN' - layout.operator("uv.cube_project") - layout.operator("uv.cylinder_project") - layout.operator("uv.sphere_project") + layout.menu("IMAGE_MT_uvs_unwrap") layout.separator() diff --git a/source/blender/blenloader/intern/versioning_300.cc b/source/blender/blenloader/intern/versioning_300.cc index cb6d2c31274..081785aa24a 100644 --- a/source/blender/blenloader/intern/versioning_300.cc +++ b/source/blender/blenloader/intern/versioning_300.cc @@ -4548,4 +4548,12 @@ void blo_do_versions_300(FileData *fd, Library * /*lib*/, Main *bmain) * * \note Keep this message at the bottom of the function. */ + { + /* Keep this block, even when empty. */ + LISTBASE_FOREACH (Scene *, scene, &bmain->scenes) { + scene->toolsettings->uvcalc_iterations = 10; + scene->toolsettings->uvcalc_weight_factor = 1.0f; + STRNCPY(scene->toolsettings->uvcalc_weight_group, "uv_importance"); + } + } } diff --git a/source/blender/editors/geometry/node_group_operator.cc b/source/blender/editors/geometry/node_group_operator.cc index b3d70eaf11c..5aded597f91 100644 --- a/source/blender/editors/geometry/node_group_operator.cc +++ b/source/blender/editors/geometry/node_group_operator.cc @@ -1142,6 +1142,7 @@ static Set get_builtin_menus(const ObjectType object_type, const eO menus.add_new("Face"); menus.add_new("Face/Face Data"); menus.add_new("UV"); + menus.add_new("UV/Unwrap"); break; case OB_MODE_SCULPT: menus.add_new("View"); diff --git a/source/blender/editors/include/ED_uvedit.hh b/source/blender/editors/include/ED_uvedit.hh index e5275135066..84897e8101c 100644 --- a/source/blender/editors/include/ED_uvedit.hh +++ b/source/blender/editors/include/ED_uvedit.hh @@ -29,6 +29,7 @@ struct bContext; struct bNode; struct bNodeTree; struct wmKeyConfig; +struct wmTimer; /* `uvedit_ops.cc` */ @@ -259,8 +260,11 @@ char ED_uvedit_select_mode_get(const Scene *scene); void ED_uvedit_select_sync_flush(const ToolSettings *ts, BMEditMesh *em, bool select); /* `uvedit_unwrap_ops.cc` */ - -void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit); +const wmTimer *ED_uvedit_live_unwrap_timer(); +/** + * \param win_modal: Support interactive (modal) unwrapping that updates with a timer. + */ +void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit, struct wmWindow *win_modal); void ED_uvedit_live_unwrap_re_solve(); void ED_uvedit_live_unwrap_end(short cancel); diff --git a/source/blender/editors/sculpt_paint/sculpt_uv.cc b/source/blender/editors/sculpt_paint/sculpt_uv.cc index 72496d516fd..9bc0334e89f 100644 --- a/source/blender/editors/sculpt_paint/sculpt_uv.cc +++ b/source/blender/editors/sculpt_paint/sculpt_uv.cc @@ -905,7 +905,8 @@ static UvSculptData *uv_sculpt_stroke_init(bContext *C, wmOperator *op, const wm data->initial_stroke->totalInitialSelected = counter; if (sima->flag & SI_LIVE_UNWRAP) { - ED_uvedit_live_unwrap_begin(scene, obedit); + wmWindow *win_modal = CTX_wm_window(C); + ED_uvedit_live_unwrap_begin(scene, obedit, win_modal); } } diff --git a/source/blender/editors/transform/transform.cc b/source/blender/editors/transform/transform.cc index 43dd24deb60..07aa2777e6d 100644 --- a/source/blender/editors/transform/transform.cc +++ b/source/blender/editors/transform/transform.cc @@ -31,6 +31,7 @@ #include "ED_node.hh" #include "ED_screen.hh" #include "ED_space_api.hh" +#include "ED_uvedit.hh" #include "ANIM_keyframing.hh" @@ -1003,6 +1004,11 @@ int transformEvent(TransInfo *t, wmOperator *op, const wmEvent *event) t->redraw |= TREDRAW_HARD; handled = true; } + else if (event->type == TIMER) { + if (ED_uvedit_live_unwrap_timer() == event->customdata) { + t->redraw |= TREDRAW_HARD; + } + } else if (!is_navigating && event->type == MOUSEMOVE) { t->mval = float2(event->mval); diff --git a/source/blender/editors/transform/transform_convert_mesh_uv.cc b/source/blender/editors/transform/transform_convert_mesh_uv.cc index 98b7dc9c98b..6e0ca4892de 100644 --- a/source/blender/editors/transform/transform_convert_mesh_uv.cc +++ b/source/blender/editors/transform/transform_convert_mesh_uv.cc @@ -390,7 +390,8 @@ static void createTransUVs(bContext *C, TransInfo *t) } if (sima->flag & SI_LIVE_UNWRAP) { - ED_uvedit_live_unwrap_begin(t->scene, tc->obedit); + wmWindow *win_modal = CTX_wm_window(C); + ED_uvedit_live_unwrap_begin(t->scene, tc->obedit, win_modal); } finally: diff --git a/source/blender/editors/uvedit/CMakeLists.txt b/source/blender/editors/uvedit/CMakeLists.txt index a9b131e2d96..87778e28a5c 100644 --- a/source/blender/editors/uvedit/CMakeLists.txt +++ b/source/blender/editors/uvedit/CMakeLists.txt @@ -44,6 +44,13 @@ set(LIB PRIVATE bf::intern::guardedalloc ) +if(WITH_UV_SLIM) + list(APPEND LIB + bf_intern_slim + ) + add_definitions(-DWITH_UV_SLIM) +endif() + blender_add_lib(bf_editor_uvedit "${SRC}" "${INC}" "${INC_SYS}" "${LIB}") diff --git a/source/blender/editors/uvedit/uvedit_ops.cc b/source/blender/editors/uvedit/uvedit_ops.cc index 7a9c9afe0e3..68f76d6be28 100644 --- a/source/blender/editors/uvedit/uvedit_ops.cc +++ b/source/blender/editors/uvedit/uvedit_ops.cc @@ -183,7 +183,7 @@ void ED_object_assign_active_image(Main *bmain, Object *ob, int mat_nr, Image *i void uvedit_live_unwrap_update(SpaceImage *sima, Scene *scene, Object *obedit) { if (sima && (sima->flag & SI_LIVE_UNWRAP)) { - ED_uvedit_live_unwrap_begin(scene, obedit); + ED_uvedit_live_unwrap_begin(scene, obedit, nullptr); ED_uvedit_live_unwrap_re_solve(); ED_uvedit_live_unwrap_end(0); } diff --git a/source/blender/editors/uvedit/uvedit_unwrap_ops.cc b/source/blender/editors/uvedit/uvedit_unwrap_ops.cc index 08438b1b5cd..6ae8f6a3686 100644 --- a/source/blender/editors/uvedit/uvedit_unwrap_ops.cc +++ b/source/blender/editors/uvedit/uvedit_unwrap_ops.cc @@ -12,8 +12,11 @@ #include "MEM_guardedalloc.h" +#include "DNA_defaults.h" +#include "DNA_meshdata_types.h" #include "DNA_modifier_types.h" #include "DNA_object_types.h" +#include "DNA_scene_defaults.h" #include "DNA_scene_types.h" #include "BKE_global.hh" @@ -36,10 +39,13 @@ #include "BKE_context.hh" #include "BKE_customdata.hh" +#include "BKE_deform.hh" #include "BKE_editmesh.hh" +#include "BKE_global.hh" #include "BKE_image.h" #include "BKE_layer.hh" #include "BKE_lib_id.hh" +#include "BKE_main.hh" #include "BKE_mesh.hh" #include "BKE_object_types.hh" #include "BKE_report.hh" @@ -74,31 +80,12 @@ using blender::Span; using blender::Vector; using blender::geometry::ParamHandle; using blender::geometry::ParamKey; +using blender::geometry::ParamSlimOptions; /* -------------------------------------------------------------------- */ /** \name Utility Functions * \{ */ -static void modifier_unwrap_state(Object *obedit, const Scene *scene, bool *r_use_subsurf) -{ - ModifierData *md; - bool subsurf = (scene->toolsettings->uvcalc_flag & UVCALC_USESUBSURF) != 0; - - md = static_cast(obedit->modifiers.first); - - /* subsurf will take the modifier settings only if modifier is first or right after mirror */ - if (subsurf) { - if (md && md->type == eModifierType_Subsurf) { - subsurf = true; - } - else { - subsurf = false; - } - } - - *r_use_subsurf = subsurf; -} - static bool ED_uvedit_ensure_uvs(Object *obedit) { if (ED_uvedit_test(obedit)) { @@ -202,6 +189,15 @@ struct UnwrapOptions { bool correct_aspect; /** Treat unselected uvs as if they were pinned. */ bool pin_unselected; + + int method; + bool use_slim; + bool use_abf; + bool use_subsurf; + bool use_weights; + + ParamSlimOptions slim; + char weight_group[MAX_VGROUP_NAME]; }; void blender::geometry::UVPackIsland_Params::setFromUnwrapOptions(const UnwrapOptions &options) @@ -213,6 +209,214 @@ void blender::geometry::UVPackIsland_Params::setFromUnwrapOptions(const UnwrapOp pin_unselected = options.pin_unselected; } +static void modifier_unwrap_state(Object *obedit, + const UnwrapOptions *options, + bool *r_use_subsurf) +{ + ModifierData *md; + bool subsurf = options->use_subsurf; + + md = static_cast(obedit->modifiers.first); + + /* Subdivision-surface will take the modifier settings + * only if modifier is first or right after mirror. */ + if (subsurf) { + if (md && md->type == eModifierType_Subsurf) { + subsurf = true; + } + else { + subsurf = false; + } + } + + *r_use_subsurf = subsurf; +} + +static UnwrapOptions unwrap_options_get(wmOperator *op, Object *ob, const ToolSettings *ts) +{ + UnwrapOptions options{}; + + /* To be set by the upper layer */ + options.topology_from_uvs = false; + options.topology_from_uvs_use_seams = false; + options.only_selected_faces = false; + options.only_selected_uvs = false; + options.pin_unselected = false; + + options.slim.skip_init = false; + + if (ts) { + options.method = ts->unwrapper; + options.correct_aspect = (ts->uvcalc_flag & UVCALC_NO_ASPECT_CORRECT) == 0; + options.fill_holes = (ts->uvcalc_flag & UVCALC_FILLHOLES) != 0; + options.use_subsurf = (ts->uvcalc_flag & UVCALC_USESUBSURF) != 0; + + options.use_weights = ts->uvcalc_flag & UVCALC_UNWRAP_USE_WEIGHTS; + STRNCPY(options.weight_group, ts->uvcalc_weight_group); + options.slim.weight_influence = ts->uvcalc_weight_factor; + + options.slim.iterations = ts->uvcalc_iterations; + options.slim.no_flip = ts->uvcalc_flag & UVCALC_UNWRAP_NO_FLIP; + } + else { + options.method = RNA_enum_get(op->ptr, "method"); + options.correct_aspect = RNA_boolean_get(op->ptr, "correct_aspect"); + options.fill_holes = RNA_boolean_get(op->ptr, "fill_holes"); + options.use_subsurf = RNA_boolean_get(op->ptr, "use_subsurf_data"); + + options.use_weights = RNA_boolean_get(op->ptr, "use_weights"); + RNA_string_get(op->ptr, "weight_group", options.weight_group); + options.slim.weight_influence = RNA_float_get(op->ptr, "weight_factor"); + + options.slim.iterations = RNA_int_get(op->ptr, "iterations"); + options.slim.no_flip = RNA_boolean_get(op->ptr, "no_flip"); + } + +#ifndef WITH_UV_SLIM + if (options.method == UVCALC_UNWRAP_METHOD_MINIMUM_STRETCH) { + options.method = UVCALC_UNWRAP_METHOD_CONFORMAL; + if (op) { + BKE_report(op->reports, RPT_WARNING, "Built without SLIM, falling back to conformal method"); + } + } +#endif /* !WITH_UV_SLIM */ + + if (options.weight_group[0] == '\0' || options.use_weights == false) { + options.slim.weight_influence = 0.0f; + } + + options.use_abf = options.method == UVCALC_UNWRAP_METHOD_ANGLE; + options.use_slim = options.method == UVCALC_UNWRAP_METHOD_MINIMUM_STRETCH; + + /* SLIM requires hole filling */ + if (options.use_slim) { + options.fill_holes = true; + } + + if (ob) { + bool use_subsurf_final; + modifier_unwrap_state(ob, &options, &use_subsurf_final); + options.use_subsurf = use_subsurf_final; + } + + return options; +} + +/* Generic sync functions + * + * NOTE: these could be moved to a generic API. + */ + +static bool rna_property_sync_flag( + PointerRNA *ptr, const char *prop_name, char flag, bool flipped, char *value_p) +{ + if (PropertyRNA *prop = RNA_struct_find_property(ptr, prop_name)) { + if (RNA_property_is_set(ptr, prop)) { + if (RNA_property_boolean_get(ptr, prop) ^ flipped) { + *value_p |= flag; + } + else { + *value_p &= ~flag; + } + return true; + } + RNA_property_boolean_set(ptr, prop, ((*value_p & flag) > 0) ^ flipped); + return false; + } + BLI_assert_unreachable(); + return false; +} + +static bool rna_property_sync_enum(PointerRNA *ptr, const char *prop_name, int *value_p) +{ + if (PropertyRNA *prop = RNA_struct_find_property(ptr, prop_name)) { + if (RNA_property_is_set(ptr, prop)) { + *value_p = RNA_property_enum_get(ptr, prop); + return true; + } + RNA_property_enum_set(ptr, prop, *value_p); + return false; + } + BLI_assert_unreachable(); + return false; +} + +static bool rna_property_sync_enum_char(PointerRNA *ptr, const char *prop_name, char *value_p) +{ + int value_i = *value_p; + if (rna_property_sync_enum(ptr, prop_name, &value_i)) { + *value_p = value_i; + return true; + } + return false; +} + +static bool rna_property_sync_int(PointerRNA *ptr, const char *prop_name, int *value_p) +{ + if (PropertyRNA *prop = RNA_struct_find_property(ptr, prop_name)) { + if (RNA_property_is_set(ptr, prop)) { + *value_p = RNA_property_int_get(ptr, prop); + return true; + } + RNA_property_int_set(ptr, prop, *value_p); + return false; + } + BLI_assert_unreachable(); + return false; +} + +static bool rna_property_sync_float(PointerRNA *ptr, const char *prop_name, float *value_p) +{ + if (PropertyRNA *prop = RNA_struct_find_property(ptr, prop_name)) { + if (RNA_property_is_set(ptr, prop)) { + *value_p = RNA_property_float_get(ptr, prop); + return true; + } + RNA_property_float_set(ptr, prop, *value_p); + return false; + } + BLI_assert_unreachable(); + return false; +} + +static bool rna_property_sync_string(PointerRNA *ptr, const char *prop_name, char value_p[]) +{ + if (PropertyRNA *prop = RNA_struct_find_property(ptr, prop_name)) { + if (RNA_property_is_set(ptr, prop)) { + RNA_property_string_get(ptr, prop, value_p); + return true; + } + RNA_property_string_set(ptr, prop, value_p); + return false; + } + BLI_assert_unreachable(); + return false; +} + +static void unwrap_options_sync_toolsettings(wmOperator *op, ToolSettings *ts) +{ + /* Remember last method for live unwrap. */ + rna_property_sync_enum_char(op->ptr, "method", &ts->unwrapper); + + /* Remember packing margin. */ + rna_property_sync_float(op->ptr, "margin", &ts->uvcalc_margin); + + rna_property_sync_int(op->ptr, "iterations", &ts->uvcalc_iterations); + + rna_property_sync_float(op->ptr, "weight_factor", &ts->uvcalc_weight_factor); + + rna_property_sync_string(op->ptr, "weight_group", ts->uvcalc_weight_group); + + rna_property_sync_flag(op->ptr, "fill_holes", UVCALC_FILLHOLES, false, &ts->uvcalc_flag); + rna_property_sync_flag( + op->ptr, "correct_aspect", UVCALC_NO_ASPECT_CORRECT, true, &ts->uvcalc_flag); + rna_property_sync_flag(op->ptr, "use_subsurf_data", UVCALC_USESUBSURF, false, &ts->uvcalc_flag); + rna_property_sync_flag(op->ptr, "no_flip", UVCALC_UNWRAP_NO_FLIP, false, &ts->uvcalc_flag); + + rna_property_sync_flag( + op->ptr, "use_weights", UVCALC_UNWRAP_USE_WEIGHTS, false, &ts->uvcalc_flag); +} + static bool uvedit_have_selection(const Scene *scene, BMEditMesh *em, const UnwrapOptions *options) { BMFace *efa; @@ -360,13 +564,17 @@ static void construct_param_handle_face_add(ParamHandle *handle, BMFace *efa, blender::geometry::ParamKey face_index, const UnwrapOptions *options, - const BMUVOffsets offsets) + const BMUVOffsets offsets, + const int cd_weight_offset, + const int cd_weight_index) { blender::Array vkeys(efa->len); blender::Array pin(efa->len); blender::Array select(efa->len); blender::Array co(efa->len); blender::Array uv(efa->len); + blender::Array weight(efa->len); + int i; BMIter liter; @@ -385,10 +593,26 @@ static void construct_param_handle_face_add(ParamHandle *handle, if (options->pin_unselected && !select[i]) { pin[i] = true; } + + /* Optional vertex group weighting. */ + if (cd_weight_offset >= 0 && cd_weight_index >= 0) { + MDeformVert *dv = (MDeformVert *)BM_ELEM_CD_GET_VOID_P(l->v, cd_weight_offset); + weight[i] = BKE_defvert_find_weight(dv, cd_weight_index); + } + else { + weight[i] = 1.0f; + } } - blender::geometry::uv_parametrizer_face_add( - handle, face_index, i, vkeys.data(), co.data(), uv.data(), pin.data(), select.data()); + blender::geometry::uv_parametrizer_face_add(handle, + face_index, + i, + vkeys.data(), + co.data(), + uv.data(), + weight.data(), + pin.data(), + select.data()); } /* Set seams on UV Parametrizer based on options. */ @@ -454,6 +678,9 @@ static ParamHandle *construct_param_handle(const Scene *scene, BM_mesh_elem_index_ensure(bm, BM_VERT); const BMUVOffsets offsets = BM_uv_map_get_offsets(bm); + const int cd_weight_offset = CustomData_get_offset(&bm->vdata, CD_MDEFORMVERT); + const int cd_weight_index = BKE_object_defgroup_name_index(ob, options->weight_group); + BM_ITER_MESH_INDEX (efa, &iter, bm, BM_FACES_OF_MESH, i) { if (uvedit_is_face_affected(scene, efa, options, offsets)) { uvedit_prepare_pinned_indices(handle, scene, efa, options, offsets); @@ -462,7 +689,8 @@ static ParamHandle *construct_param_handle(const Scene *scene, BM_ITER_MESH_INDEX (efa, &iter, bm, BM_FACES_OF_MESH, i) { if (uvedit_is_face_affected(scene, efa, options, offsets)) { - construct_param_handle_face_add(handle, scene, efa, i, options, offsets); + construct_param_handle_face_add( + handle, scene, efa, i, options, offsets, cd_weight_offset, cd_weight_index); } } @@ -507,6 +735,9 @@ static ParamHandle *construct_param_handle_multi(const Scene *scene, continue; } + const int cd_weight_offset = CustomData_get_offset(&bm->vdata, CD_MDEFORMVERT); + const int cd_weight_index = BKE_object_defgroup_name_index(obedit, options->weight_group); + BM_ITER_MESH_INDEX (efa, &iter, bm, BM_FACES_OF_MESH, i) { if (uvedit_is_face_affected(scene, efa, options, offsets)) { uvedit_prepare_pinned_indices(handle, scene, efa, options, offsets); @@ -515,7 +746,8 @@ static ParamHandle *construct_param_handle_multi(const Scene *scene, BM_ITER_MESH_INDEX (efa, &iter, bm, BM_FACES_OF_MESH, i) { if (uvedit_is_face_affected(scene, efa, options, offsets)) { - construct_param_handle_face_add(handle, scene, efa, i + offset, options, offsets); + construct_param_handle_face_add( + handle, scene, efa, i + offset, options, offsets, cd_weight_offset, cd_weight_index); } } @@ -608,6 +840,7 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene, BMEdge **edgeMap; const BMUVOffsets offsets = BM_uv_map_get_offsets(em->bm); + const int cd_weight_index = BKE_object_defgroup_name_index(ob, options->weight_group); ParamHandle *handle = new blender::geometry::ParamHandle(); @@ -630,6 +863,7 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene, const blender::Span subsurf_edges = subdiv_mesh->edges(); const blender::OffsetIndices subsurf_facess = subdiv_mesh->faces(); const blender::Span subsurf_corner_verts = subdiv_mesh->corner_verts(); + const blender::Span subsurf_deform_verts = subdiv_mesh->deform_verts(); const int *origVertIndices = static_cast( CustomData_get_layer(&subdiv_mesh->vert_data, CD_ORIGINDEX)); @@ -666,6 +900,7 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene, bool pin[4], select[4]; const float *co[4]; float *uv[4]; + float weight[4]; BMFace *origFace = faceMap[i]; if (scene->toolsettings->uv_flag & UV_SYNC_SELECTION) { @@ -696,6 +931,24 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene, co[2] = subsurf_positions[poly_corner_verts[2]]; co[3] = subsurf_positions[poly_corner_verts[3]]; + /* Optional vertex group weights. */ + if (cd_weight_index >= 0) { + weight[0] = BKE_defvert_find_weight(&subsurf_deform_verts[poly_corner_verts[0]], + cd_weight_index); + weight[1] = BKE_defvert_find_weight(&subsurf_deform_verts[poly_corner_verts[1]], + cd_weight_index); + weight[2] = BKE_defvert_find_weight(&subsurf_deform_verts[poly_corner_verts[2]], + cd_weight_index); + weight[3] = BKE_defvert_find_weight(&subsurf_deform_verts[poly_corner_verts[3]], + cd_weight_index); + } + else { + weight[0] = 1.0f; + weight[1] = 1.0f; + weight[2] = 1.0f; + weight[3] = 1.0f; + } + /* This is where all the magic is done. * If the vertex exists in the, we pass the original uv pointer to the solver, thus * flushing the solution to the edit mesh. */ @@ -728,7 +981,8 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene, &pin[3], &select[3]); - blender::geometry::uv_parametrizer_face_add(handle, key, 4, vkeys, co, uv, pin, select); + blender::geometry::uv_parametrizer_face_add( + handle, key, 4, vkeys, co, uv, weight, pin, select); } /* These are calculated from original mesh too. */ @@ -1434,7 +1688,7 @@ static int pack_islands_exec(bContext *C, wmOperator *op) const Scene *scene = CTX_data_scene(C); const SpaceImage *sima = CTX_wm_space_image(C); - UnwrapOptions options{}; + UnwrapOptions options = unwrap_options_get(op, nullptr, nullptr); options.topology_from_uvs = true; options.only_selected_faces = true; options.only_selected_uvs = true; @@ -1725,7 +1979,7 @@ static int average_islands_scale_exec(bContext *C, wmOperator *op) ToolSettings *ts = scene->toolsettings; const bool synced_selection = (ts->uv_flag & UV_SYNC_SELECTION) != 0; - UnwrapOptions options{}; + UnwrapOptions options = unwrap_options_get(nullptr, nullptr, ts); options.topology_from_uvs = true; options.only_selected_faces = true; options.only_selected_uvs = true; @@ -1788,36 +2042,49 @@ void UV_OT_average_islands_scale(wmOperatorType *ot) static struct { ParamHandle **handles; uint len, len_alloc; + wmTimer *timer; } g_live_unwrap = {nullptr}; -void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit) +const wmTimer *ED_uvedit_live_unwrap_timer() +{ + return g_live_unwrap.timer; +} + +void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit, wmWindow *win_modal) { ParamHandle *handle = nullptr; BMEditMesh *em = BKE_editmesh_from_object(obedit); - const bool abf = (scene->toolsettings->unwrapper == UVCALC_UNWRAP_METHOD_ANGLE); - bool use_subsurf; - - modifier_unwrap_state(obedit, scene, &use_subsurf); if (!ED_uvedit_test(obedit)) { return; } - UnwrapOptions options{}; + UnwrapOptions options = unwrap_options_get(nullptr, obedit, scene->toolsettings); options.topology_from_uvs = false; options.only_selected_faces = false; options.only_selected_uvs = false; - options.fill_holes = (scene->toolsettings->uvcalc_flag & UVCALC_FILLHOLES) != 0; - options.correct_aspect = (scene->toolsettings->uvcalc_flag & UVCALC_NO_ASPECT_CORRECT) == 0; - if (use_subsurf) { + if (options.use_subsurf) { handle = construct_param_handle_subsurfed(scene, obedit, em, &options, nullptr); } else { handle = construct_param_handle(scene, obedit, em->bm, &options, nullptr); } - blender::geometry::uv_parametrizer_lscm_begin(handle, true, abf); + if (options.use_slim) { + options.slim.no_flip = false; + options.slim.skip_init = true; + uv_parametrizer_slim_live_begin(handle, &options.slim); + + if (win_modal) { + wmWindowManager *wm = static_cast(G_MAIN->wm.first); + BLI_assert(!g_live_unwrap.timer); + g_live_unwrap.timer = WM_event_timer_add(wm, win_modal, TIMER, 0.01f); + } + } + else { + blender::geometry::uv_parametrizer_lscm_begin(handle, true, options.use_abf); + } /* Create or increase size of g_live_unwrap.handles array */ if (g_live_unwrap.handles == nullptr) { @@ -1839,7 +2106,13 @@ void ED_uvedit_live_unwrap_re_solve() { if (g_live_unwrap.handles) { for (int i = 0; i < g_live_unwrap.len; i++) { - blender::geometry::uv_parametrizer_lscm_solve(g_live_unwrap.handles[i], nullptr, nullptr); + if (uv_parametrizer_is_slim(g_live_unwrap.handles[i])) { + uv_parametrizer_slim_live_solve_iteration(g_live_unwrap.handles[i]); + } + else { + blender::geometry::uv_parametrizer_lscm_solve(g_live_unwrap.handles[i], nullptr, nullptr); + } + blender::geometry::uv_parametrizer_flush(g_live_unwrap.handles[i]); } } @@ -1847,9 +2120,22 @@ void ED_uvedit_live_unwrap_re_solve() void ED_uvedit_live_unwrap_end(short cancel) { + if (g_live_unwrap.timer) { + wmWindow *win = g_live_unwrap.timer->win; + wmWindowManager *wm = static_cast(G_MAIN->wm.first); + WM_event_timer_remove(wm, win, g_live_unwrap.timer); + g_live_unwrap.timer = nullptr; + } + if (g_live_unwrap.handles) { for (int i = 0; i < g_live_unwrap.len; i++) { - blender::geometry::uv_parametrizer_lscm_end(g_live_unwrap.handles[i]); + if (uv_parametrizer_is_slim(g_live_unwrap.handles[i])) { + uv_parametrizer_slim_live_end(g_live_unwrap.handles[i]); + } + else { + blender::geometry::uv_parametrizer_lscm_end(g_live_unwrap.handles[i]); + } + if (cancel) { blender::geometry::uv_parametrizer_flush_restore(g_live_unwrap.handles[i]); } @@ -2365,7 +2651,7 @@ static void uvedit_unwrap(const Scene *scene, } bool use_subsurf; - modifier_unwrap_state(obedit, scene, &use_subsurf); + modifier_unwrap_state(obedit, options, &use_subsurf); ParamHandle *handle; if (use_subsurf) { @@ -2375,10 +2661,14 @@ static void uvedit_unwrap(const Scene *scene, handle = construct_param_handle(scene, obedit, em->bm, options, r_count_failed); } - blender::geometry::uv_parametrizer_lscm_begin( - handle, false, scene->toolsettings->unwrapper == UVCALC_UNWRAP_METHOD_ANGLE); - blender::geometry::uv_parametrizer_lscm_solve(handle, r_count_changed, r_count_failed); - blender::geometry::uv_parametrizer_lscm_end(handle); + if (options->use_slim) { + uv_parametrizer_slim_solve(handle, &options->slim, r_count_changed, r_count_failed); + } + else { + blender::geometry::uv_parametrizer_lscm_begin(handle, false, options->use_abf); + blender::geometry::uv_parametrizer_lscm_solve(handle, r_count_changed, r_count_failed); + blender::geometry::uv_parametrizer_lscm_end(handle); + } blender::geometry::uv_parametrizer_average(handle, true, false, false); @@ -2403,7 +2693,7 @@ static void uvedit_unwrap_multi(const Scene *scene, void ED_uvedit_live_unwrap(const Scene *scene, const Span objects) { if (scene->toolsettings->edge_mode_live_unwrap) { - UnwrapOptions options{}; + UnwrapOptions options = unwrap_options_get(nullptr, nullptr, scene->toolsettings); options.topology_from_uvs = false; options.only_selected_faces = false; options.only_selected_uvs = false; @@ -2431,23 +2721,24 @@ static int unwrap_exec(bContext *C, wmOperator *op) { ViewLayer *view_layer = CTX_data_view_layer(C); const Scene *scene = CTX_data_scene(C); - int method = RNA_enum_get(op->ptr, "method"); - const bool use_subsurf = RNA_boolean_get(op->ptr, "use_subsurf_data"); int reported_errors = 0; - /* We will report an error unless at least one object - * has the subsurf modifier in the right place. */ - bool subsurf_error = use_subsurf; Vector objects = BKE_view_layer_array_from_objects_in_edit_mode_unique_data( scene, view_layer, CTX_wm_view3d(C)); - UnwrapOptions options{}; + unwrap_options_sync_toolsettings(op, scene->toolsettings); + + UnwrapOptions options = unwrap_options_get(op, nullptr, nullptr); options.topology_from_uvs = false; options.only_selected_faces = true; options.only_selected_uvs = false; options.fill_holes = RNA_boolean_get(op->ptr, "fill_holes"); options.correct_aspect = RNA_boolean_get(op->ptr, "correct_aspect"); + /* We will report an error unless at least one object + * has the subsurf modifier in the right place. */ + bool subsurf_error = options.use_subsurf; + if (CTX_wm_space_image(C)) { /* Inside the UV Editor, only unwrap selected UVs. */ options.only_selected_uvs = true; @@ -2470,7 +2761,7 @@ static int unwrap_exec(bContext *C, wmOperator *op) if (subsurf_error) { /* Double up the check here but better keep uvedit_unwrap interface simple and not * pass operator for warning append. */ - modifier_unwrap_state(obedit, scene, &use_subsurf_final); + modifier_unwrap_state(obedit, &options, &use_subsurf_final); if (use_subsurf_final) { subsurf_error = false; } @@ -2507,43 +2798,6 @@ static int unwrap_exec(bContext *C, wmOperator *op) "Subdivision Surface modifier needs to be first to work with unwrap"); } - /* remember last method for live unwrap */ - if (RNA_struct_property_is_set(op->ptr, "method")) { - scene->toolsettings->unwrapper = method; - } - else { - RNA_enum_set(op->ptr, "method", scene->toolsettings->unwrapper); - } - - /* remember packing margin */ - if (RNA_struct_property_is_set(op->ptr, "margin")) { - scene->toolsettings->uvcalc_margin = RNA_float_get(op->ptr, "margin"); - } - else { - RNA_float_set(op->ptr, "margin", scene->toolsettings->uvcalc_margin); - } - - if (options.fill_holes) { - scene->toolsettings->uvcalc_flag |= UVCALC_FILLHOLES; - } - else { - scene->toolsettings->uvcalc_flag &= ~UVCALC_FILLHOLES; - } - - if (options.correct_aspect) { - scene->toolsettings->uvcalc_flag &= ~UVCALC_NO_ASPECT_CORRECT; - } - else { - scene->toolsettings->uvcalc_flag |= UVCALC_NO_ASPECT_CORRECT; - } - - if (use_subsurf) { - scene->toolsettings->uvcalc_flag |= UVCALC_USESUBSURF; - } - else { - scene->toolsettings->uvcalc_flag &= ~UVCALC_USESUBSURF; - } - /* execute unwrap */ int count_changed = 0; int count_failed = 0; @@ -2575,11 +2829,56 @@ static int unwrap_exec(bContext *C, wmOperator *op) return OPERATOR_FINISHED; } +static void unwrap_draw(bContext * /*C*/, wmOperator *op) +{ + uiLayout *layout = op->layout; + + uiLayoutSetPropSep(layout, true); + uiLayoutSetPropDecorate(layout, false); + + /* Main draw call */ + PointerRNA ptr = RNA_pointer_create(nullptr, op->type->srna, op->properties); + + uiLayout *col; + + col = uiLayoutColumn(layout, true); + uiItemR(col, &ptr, "method", UI_ITEM_NONE, nullptr, ICON_NONE); + bool is_slim = RNA_enum_get(op->ptr, "method") == UVCALC_UNWRAP_METHOD_MINIMUM_STRETCH; + + if (is_slim) { + uiItemR(col, &ptr, "iterations", UI_ITEM_NONE, nullptr, ICON_NONE); + uiItemR(col, &ptr, "no_flip", UI_ITEM_NONE, nullptr, ICON_NONE); + + uiItemS(col); + uiItemR(col, &ptr, "use_weights", UI_ITEM_NONE, nullptr, ICON_NONE); + + if (RNA_boolean_get(op->ptr, "use_weights")) { + col = uiLayoutColumn(layout, true); + uiItemR(col, &ptr, "weight_group", UI_ITEM_NONE, nullptr, ICON_NONE); + uiItemR(col, &ptr, "weight_factor", UI_ITEM_NONE, nullptr, ICON_NONE); + } + } + else { + uiItemR(col, &ptr, "fill_holes", UI_ITEM_NONE, nullptr, ICON_NONE); + } + + uiItemS(col); + uiItemR(col, &ptr, "use_subsurf_data", UI_ITEM_NONE, nullptr, ICON_NONE); + + uiItemS(col); + uiItemR(col, &ptr, "correct_aspect", UI_ITEM_NONE, nullptr, ICON_NONE); + uiItemR(col, &ptr, "margin_method", UI_ITEM_NONE, nullptr, ICON_NONE); + uiItemR(col, &ptr, "margin", UI_ITEM_NONE, nullptr, ICON_NONE); +} + void UV_OT_unwrap(wmOperatorType *ot) { + const ToolSettings *tool_settings_default = DNA_struct_default_get(ToolSettings); + static const EnumPropertyItem method_items[] = { {UVCALC_UNWRAP_METHOD_ANGLE, "ANGLE_BASED", 0, "Angle Based", ""}, {UVCALC_UNWRAP_METHOD_CONFORMAL, "CONFORMAL", 0, "Conformal", ""}, + {UVCALC_UNWRAP_METHOD_MINIMUM_STRETCH, "MINIMUM_STRETCH", 0, "Minimum Stretch", ""}, {0, nullptr, 0, nullptr, nullptr}, }; @@ -2593,23 +2892,27 @@ void UV_OT_unwrap(wmOperatorType *ot) ot->exec = unwrap_exec; ot->poll = ED_operator_uvmap; + /* Only draw relevant ui elements */ + ot->ui = unwrap_draw; + /* properties */ - RNA_def_enum(ot->srna, - "method", - method_items, - 0, - "Method", - "Unwrapping method (Angle Based usually gives better results than Conformal, while " - "being somewhat slower)"); + ot->prop = RNA_def_enum( + ot->srna, + "method", + method_items, + tool_settings_default->unwrapper, + "Method", + "Unwrapping method (Angle Based usually gives better results than Conformal, while " + "being somewhat slower)"); RNA_def_boolean(ot->srna, "fill_holes", - true, + tool_settings_default->uvcalc_flag & UVCALC_FILLHOLES, "Fill Holes", "Virtually fill holes in mesh before unwrapping, to better avoid overlaps and " "preserve symmetry"); RNA_def_boolean(ot->srna, "correct_aspect", - true, + !(tool_settings_default->uvcalc_flag & UVCALC_NO_ASPECT_CORRECT), "Correct Aspect", "Map UVs taking image aspect ratio into account"); RNA_def_boolean( @@ -2626,6 +2929,46 @@ void UV_OT_unwrap(wmOperatorType *ot) ""); RNA_def_float_factor( ot->srna, "margin", 0.001f, 0.0f, 1.0f, "Margin", "Space between islands", 0.0f, 1.0f); + + /* SLIM only */ + RNA_def_boolean(ot->srna, + "no_flip", + tool_settings_default->uvcalc_flag & UVCALC_UNWRAP_NO_FLIP, + "No Flip", + "Prevent flipping UV's, " + "flipping may lower distortion depending on the position of pins"); + + RNA_def_int(ot->srna, + "iterations", + tool_settings_default->uvcalc_iterations, + 0, + 10000, + "Iterations", + "Number of iterations when \"Minimum Stretch\" method is used", + 1, + 30); + + RNA_def_boolean(ot->srna, + "use_weights", + tool_settings_default->uvcalc_flag & UVCALC_UNWRAP_USE_WEIGHTS, + "Importance Weights", + "Whether to take into account per-vertex importance weights"); + RNA_def_string(ot->srna, + "weight_group", + tool_settings_default->uvcalc_weight_group, + MAX_ID_NAME, + "Weight Group", + "Vertex group name for importance weights (modulating the deform)"); + RNA_def_float( + ot->srna, + "weight_factor", + tool_settings_default->uvcalc_weight_factor, + -10000.0, + 10000.0, + "Weight Factor", + "How much influence the weightmap has for weighted parameterization, 0 being no influence", + -10.0, + 10.0); } /** \} */ diff --git a/source/blender/geometry/CMakeLists.txt b/source/blender/geometry/CMakeLists.txt index f692c3386db..22c9aaa13d6 100644 --- a/source/blender/geometry/CMakeLists.txt +++ b/source/blender/geometry/CMakeLists.txt @@ -101,6 +101,13 @@ set(LIB PRIVATE bf::extern::fmtlib ) +if(WITH_UV_SLIM) + list(APPEND INC + ../../../intern/slim + ) + add_definitions(-DWITH_UV_SLIM) +endif() + if(WITH_OPENVDB) list(APPEND INC ../../../intern/openvdb diff --git a/source/blender/geometry/GEO_uv_parametrizer.hh b/source/blender/geometry/GEO_uv_parametrizer.hh index cd1e1b5d22b..14ac0a9173d 100644 --- a/source/blender/geometry/GEO_uv_parametrizer.hh +++ b/source/blender/geometry/GEO_uv_parametrizer.hh @@ -6,6 +6,10 @@ #include "BLI_sys_types.h" /* for intptr_t support */ +namespace slim { +struct MatrixTransfer; +} + /** \file * \ingroup geo */ @@ -55,6 +59,9 @@ class ParamHandle { RNG *rng = nullptr; float blend = 0.0f; + + /* SLIM uv unwrapping */ + slim::MatrixTransfer *slim_mt = nullptr; }; /* -------------------------------------------------------------------- */ @@ -84,6 +91,7 @@ void uv_parametrizer_face_add(ParamHandle *handle, const ParamKey *vkeys, const float **co, float **uv, /* Output will eventually be written to `uv`. */ + float *weight, const bool *pin, const bool *select); @@ -96,6 +104,34 @@ void uv_parametrizer_construct_end(ParamHandle *handle, /** \} */ +/* -------------------------------------------------------------------- */ +/** \name SLIM: + * + * - begin: data is gathered into matrices and transferred to SLIM. + * - solve: compute cheap initialization (if necessary) and refine iteratively. + * - end: clean up. + * \{ */ + +struct ParamSlimOptions { + float weight_influence = 0.0f; + int iterations = 0; + bool no_flip = false; + bool skip_init = false; +}; + +void uv_parametrizer_slim_solve(ParamHandle *phandle, + const ParamSlimOptions *slim_options, + int *count_changed, + int *count_failed); + +void uv_parametrizer_slim_live_begin(ParamHandle *phandle, const ParamSlimOptions *slim_options); +void uv_parametrizer_slim_live_solve_iteration(ParamHandle *phandle); +void uv_parametrizer_slim_live_end(ParamHandle *phandle); +void uv_parametrizer_slim_stretch_iteration(ParamHandle *phandle, float blend); +bool uv_parametrizer_is_slim(const ParamHandle *phandle); + +/** \} */ + /* -------------------------------------------------------------------- */ /** \name Least Squares Conformal Maps: * diff --git a/source/blender/geometry/intern/uv_parametrizer.cc b/source/blender/geometry/intern/uv_parametrizer.cc index f6b166f9ef0..b13123a8788 100644 --- a/source/blender/geometry/intern/uv_parametrizer.cc +++ b/source/blender/geometry/intern/uv_parametrizer.cc @@ -6,6 +6,9 @@ * \ingroup eduv */ +#include +#include + #include "GEO_uv_parametrizer.hh" #include "BLI_array.hh" @@ -19,6 +22,10 @@ #include "BLI_polyfill_2d_beautify.h" #include "BLI_rand.h" +#ifdef WITH_UV_SLIM +# include "slim_matrix_transfer.h" +#endif + #include "GEO_uv_pack.hh" #include "eigen_capi.h" @@ -30,6 +37,13 @@ namespace blender::geometry { #define param_warning(message) \ {/* `printf("Warning %s:%d: %s\n", __FILE__, __LINE__, message);` */}(void)0 +/* Prevent unused function warnings when slim is disabled. */ +#ifdef WITH_UV_SLIM +# define UNUSED_FUNCTION_NO_SLIM(x) x +#else +# define UNUSED_FUNCTION_NO_SLIM UNUSED_FUNCTION +#endif + /* Special Purpose Hash */ using PHashKey = uintptr_t; @@ -60,6 +74,10 @@ struct PVert { float co[3]; float uv[2]; uint flag; + + float weight; + bool on_boundary_flag; + int slim_id; }; struct PEdge { @@ -121,6 +139,7 @@ enum PFaceFlag { PFACE_CONNECTED = 1, PFACE_FILLED = 2, PFACE_COLLAPSE = 4, + PFACE_DONE = 8, }; /* Chart */ @@ -344,6 +363,31 @@ static void p_face_angles(PFace *f, double *r_a1, double *r_a2, double *r_a3) p_triangle_angles(v1->co, v2->co, v3->co, r_a1, r_a2, r_a3); } +static float p_vec_cos(const float v1[3], const float v2[3], const float v3[3]) +{ + return cos_v3v3v3(v1, v2, v3); +} + +static void p_triangle_cos(const float v1[3], + const float v2[3], + const float v3[3], + float *r_cos1, + float *r_cos2, + float *r_cos3) +{ + *r_cos1 = p_vec_cos(v3, v1, v2); + *r_cos2 = p_vec_cos(v1, v2, v3); + *r_cos3 = p_vec_cos(v2, v3, v1); +} + +static void UNUSED_FUNCTION(p_face_cos)(PFace *f, float *r_cos1, float *r_cos2, float *r_cos3) +{ + PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next; + PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert; + + p_triangle_cos(v1->co, v2->co, v3->co, r_cos1, r_cos2, r_cos3); +} + static float p_face_area(PFace *f) { PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next; @@ -371,6 +415,11 @@ static float p_edge_length(PEdge *e) return len_v3v3(e->vert->co, e->next->vert->co); } +static float p_edge_length_squared(PEdge *e) +{ + return len_squared_v3v3(e->vert->co, e->next->vert->co); +} + static float p_edge_uv_length(PEdge *e) { return len_v2v2(e->vert->uv, e->next->vert->uv); @@ -478,6 +527,11 @@ static PEdge *p_wheel_edge_next(PEdge *e) return e->next->next->pair; } +static const PEdge *p_wheel_edge_next(const PEdge *e) +{ + return e->next->next->pair; +} + static PEdge *p_wheel_edge_prev(PEdge *e) { return (e->pair) ? e->pair->next : nullptr; @@ -592,6 +646,8 @@ static void p_vert_load_pin_select_uvs(ParamHandle *handle, PVert *v) } } +static void p_chart_flush_collapsed_uvs(PChart *chart); + static void p_flush_uvs(ParamHandle *handle, PChart *chart) { const float blend = handle->blend; @@ -603,6 +659,17 @@ static void p_flush_uvs(ParamHandle *handle, PChart *chart) e->orig_uv[1] = blend * e->old_uv[1] + invblend * e->vert->uv[1]; } } + + if (chart->collapsed_edges) { + p_chart_flush_collapsed_uvs(chart); + + for (PEdge *e = chart->collapsed_edges; e; e = e->nextlink) { + if (e->orig_uv) { + e->orig_uv[0] = blend * e->old_uv[0] + invblend_x * e->vert->uv[0]; + e->orig_uv[1] = blend * e->old_uv[1] + invblend * e->vert->uv[1]; + } + } + } } static void p_face_backup_uvs(PFace *f) @@ -643,10 +710,12 @@ static void p_face_restore_uvs(PFace *f) /* Construction (use only during construction, relies on u.key being set */ -static PVert *p_vert_add(ParamHandle *handle, PHashKey key, const float co[3], PEdge *e) +static PVert *p_vert_add( + ParamHandle *handle, PHashKey key, const float co[3], const float weight, PEdge *e) { PVert *v = (PVert *)BLI_memarena_alloc(handle->arena, sizeof(*v)); copy_v3_v3(v->co, co); + v->weight = weight; /* Sanity check, a single nan/inf point causes the entire result to be invalid. * Note that values within the calculation may _become_ non-finite, @@ -666,14 +735,15 @@ static PVert *p_vert_add(ParamHandle *handle, PHashKey key, const float co[3], P return v; } -static PVert *p_vert_lookup(ParamHandle *handle, PHashKey key, const float co[3], PEdge *e) +static PVert *p_vert_lookup( + ParamHandle *handle, PHashKey key, const float co[3], const float weight, PEdge *e) { PVert *v = (PVert *)phash_lookup(handle->hash_verts, key); if (v) { return v; } - return p_vert_add(handle, key, co, e); + return p_vert_add(handle, key, co, weight, e); } static PVert *p_vert_copy(ParamHandle *handle, PVert *v) @@ -1014,6 +1084,7 @@ static PFace *p_face_add_construct(ParamHandle *handle, const ParamKey *vkeys, const float **co, float **uv, + const float *weight, int i1, int i2, int i3, @@ -1023,9 +1094,21 @@ static PFace *p_face_add_construct(ParamHandle *handle, PFace *f = p_face_add(handle); PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next; - e1->vert = p_vert_lookup(handle, vkeys[i1], co[i1], e1); - e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], e2); - e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], e3); + float weight1, weight2, weight3; + if (weight) { + weight1 = weight[i1]; + weight2 = weight[i2]; + weight3 = weight[i3]; + } + else { + weight1 = 1.0f; + weight2 = 1.0f; + weight3 = 1.0f; + } + + e1->vert = p_vert_lookup(handle, vkeys[i1], co[i1], weight1, e1); + e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], weight2, e2); + e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], weight3, e3); e1->orig_uv = uv[i1]; e2->orig_uv = uv[i2]; @@ -1556,6 +1639,7 @@ static void p_vert_harmonic_insert(PVert *v) p_vert_map_harmonic_weights(v); } +#endif static void p_vert_fix_edge_pointer(PVert *v) { @@ -1680,6 +1764,7 @@ static void p_collapse_edge(PEdge *edge, PEdge *pair) } } +#if 0 static void p_split_vertex(PEdge *edge, PEdge *pair) { PVert *newv, *keepv; @@ -1727,6 +1812,7 @@ static void p_split_vertex(PEdge *edge, PEdge *pair) e = p_wheel_edge_next(e); } while (e && (e != newv->edge)); } +#endif static bool p_collapse_allowed_topologic(PEdge *edge, PEdge *pair) { @@ -1753,6 +1839,7 @@ static bool p_collapse_allowed_topologic(PEdge *edge, PEdge *pair) return true; } +#if 0 static bool p_collapse_normal_flipped(float *v1, float *v2, float *vold, float *vnew) { float nold[3], nnew[3], sub1[3], sub2[3]; @@ -1944,6 +2031,7 @@ static float p_collapse_cost(PEdge *edge, PEdge *pair) return cost; } +#endif static void p_collapse_cost_vertex( PVert *vert, @@ -2066,38 +2154,6 @@ static void p_chart_simplify_compute(PChart *chart, /* For debugging. */ static const int MAX_SIMPLIFY = INT_MAX; - PVert *v, *nextv = nullptr; - PEdge *e, *nexte = nullptr; - PFace *f, *nextf = nullptr; - - for (v = chart->collapsed_verts; v; v = nextv) { - nextv = v->nextlink; - v->nextlink = chart->verts; - chart->verts = v; - chart->nverts++; - } - - for (e = chart->collapsed_edges; e; e = nexte) { - nexte = e->nextlink; - e->nextlink = chart->edges; - chart->edges = e; - chart->nedges++; - } - - for (f = chart->collapsed_faces; f; f = nextf) { - nextf = f->nextlink; - f->nextlink = chart->faces; - chart->faces = f; - chart->nfaces++; - } - - chart->collapsed_verts = nullptr; - chart->collapsed_edges = nullptr; - chart->collapsed_faces = nullptr; -} - -static void p_chart_simplify_compute(PChart *chart) -{ /* Computes a list of edge collapses / vertex splits. The collapsed * simplices go in the `chart->collapsed_*` lists, The original and * collapsed may then be view as stacks, where the next collapse/split @@ -2198,6 +2254,41 @@ static void p_chart_simplify_compute(PChart *chart) p_chart_post_collapse_flush(chart, collapsededges); } +#if 0 +static void p_chart_post_split_flush(PChart *chart) +{ + /* Move from `collapsed_*`. */ + + PVert *v, *nextv = nullptr; + PEdge *e, *nexte = nullptr; + PFace *f, *nextf = nullptr; + + for (v = chart->collapsed_verts; v; v = nextv) { + nextv = v->nextlink; + v->nextlink = chart->verts; + chart->verts = v; + chart->nverts++; + } + + for (e = chart->collapsed_edges; e; e = nexte) { + nexte = e->nextlink; + e->nextlink = chart->edges; + chart->edges = e; + chart->nedges++; + } + + for (f = chart->collapsed_faces; f; f = nextf) { + nextf = f->nextlink; + f->nextlink = chart->faces; + chart->faces = f; + chart->nfaces++; + } + + chart->collapsed_verts = nullptr; + chart->collapsed_edges = nullptr; + chart->collapsed_faces = nullptr; +} + static void p_chart_complexify(PChart *chart) { /* For debugging. */ @@ -2235,12 +2326,10 @@ static void p_chart_complexify(PChart *chart) p_chart_post_split_flush(chart); } -# if 0 static void p_chart_simplify(PChart *chart) { /* Not implemented, needs proper reordering in split_flush. */ } -# endif #endif /* ABF */ @@ -3752,6 +3841,7 @@ static void p_add_ngon(ParamHandle *handle, const ParamKey *vkeys, const float **co, float **uv, /* Output will eventually be written to `uv`. */ + const float *weight, const bool *pin, const bool *select) { @@ -3800,10 +3890,22 @@ static void p_add_ngon(ParamHandle *handle, const ParamKey tri_vkeys[3] = {vkeys[v0], vkeys[v1], vkeys[v2]}; const float *tri_co[3] = {co[v0], co[v1], co[v2]}; float *tri_uv[3] = {uv[v0], uv[v1], uv[v2]}; + + float *tri_weight = nullptr; + float tri_weight_tmp[3]; + + if (weight) { + tri_weight_tmp[0] = weight[v0]; + tri_weight_tmp[1] = weight[v1]; + tri_weight_tmp[2] = weight[v2]; + tri_weight = tri_weight_tmp; + }; + const bool tri_pin[3] = {pin[v0], pin[v1], pin[v2]}; const bool tri_select[3] = {select[v0], select[v1], select[v2]}; - uv_parametrizer_face_add(handle, key, 3, tri_vkeys, tri_co, tri_uv, tri_pin, tri_select); + uv_parametrizer_face_add( + handle, key, 3, tri_vkeys, tri_co, tri_uv, tri_weight, tri_pin, tri_select); } BLI_memarena_clear(arena); @@ -3815,6 +3917,7 @@ void uv_parametrizer_face_add(ParamHandle *phandle, const ParamKey *vkeys, const float **co, float **uv, + float *weight, const bool *pin, const bool *select) { @@ -3852,7 +3955,7 @@ void uv_parametrizer_face_add(ParamHandle *phandle, * NOTE: Should probably call `uv_parametrizer_face_add` * instead of `p_face_add_construct`. */ int iprev = permute[(i + pm - 1) % pm]; - p_face_add_construct(phandle, key, vkeys, co, uv, iprev, i0, i1, pin, select); + p_face_add_construct(phandle, key, vkeys, co, uv, weight, iprev, i0, i1, pin, select); permute.remove(i); if (permute.size() == 3) { @@ -3865,6 +3968,7 @@ void uv_parametrizer_face_add(ParamHandle *phandle, Array vkeys_sub(pm); Array co_sub(pm); Array uv_sub(pm); + Array weight_sub(pm); Array pin_sub(pm); Array select_sub(pm); for (int i = 0; i < pm; i++) { @@ -3872,6 +3976,7 @@ void uv_parametrizer_face_add(ParamHandle *phandle, vkeys_sub[i] = vkeys[j]; co_sub[i] = co[j]; uv_sub[i] = uv[j]; + weight_sub[i] = weight[j]; pin_sub[i] = pin && pin[j]; select_sub[i] = select && select[j]; } @@ -3881,6 +3986,7 @@ void uv_parametrizer_face_add(ParamHandle *phandle, &vkeys_sub.first(), &co_sub.first(), &uv_sub.first(), + &weight_sub.first(), &pin_sub.first(), &select_sub.first()); return; /* Nothing more to do. */ @@ -3889,11 +3995,11 @@ void uv_parametrizer_face_add(ParamHandle *phandle, } if (nverts > 3) { /* ngon */ - p_add_ngon(phandle, key, nverts, vkeys, co, uv, pin, select); + p_add_ngon(phandle, key, nverts, vkeys, co, uv, weight, pin, select); } else if (!p_face_exists(phandle, vkeys, 0, 1, 2)) { /* triangle */ - p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select); + p_face_add_construct(phandle, key, vkeys, co, uv, weight, 0, 1, 2, pin, select); } } @@ -4283,4 +4389,1005 @@ void uv_parametrizer_flush_restore(ParamHandle *phandle) } } +/* -------------------------------------------------------------------- */ +/** \name Degenerate Geometry Fixing + * \{ */ + +static bool p_collapse_doubles_allowed(PEdge *edge, PEdge *pair, float threshold_squared) +{ + PVert *oldv, *keepv; + + p_collapsing_verts(edge, pair, &oldv, &keepv); + + /* Do not collapse a pinned vertex unless the target vertex + * is also pinned. */ + if ((oldv->flag & PVERT_PIN) && !(keepv->flag & PVERT_PIN)) { + return false; + } + + if (!p_collapse_allowed_topologic(edge, pair)) { + return false; + } + + PEdge *collapse_e = edge ? edge : pair; + return p_edge_length_squared(collapse_e) < threshold_squared; +} + +static float p_collapse_doubles_cost(PEdge *edge, PEdge *pair) +{ + PEdge *collapse_e = edge ? edge : pair; + return p_edge_length_squared(collapse_e); +} + +static void UNUSED_FUNCTION(p_chart_collapse_doubles)(PChart *chart, const float threshold) +{ + const float threshold_squared = threshold * threshold; + + p_chart_simplify_compute(chart, p_collapse_doubles_cost, [=](PEdge *e, PEdge *pair) { + return p_collapse_doubles_allowed(e, pair, threshold_squared); + }); +} + +static void p_chart_flush_collapsed_uvs(PChart *chart) +{ + PEdge *e, *pair, *edge; + PVert *newv, *keepv; + + for (e = chart->collapsed_edges; e; e = e->nextlink) { + if (!(e->flag & PEDGE_COLLAPSE_EDGE)) { + break; + } + edge = e; + pair = e->pair; + if (edge->flag & PEDGE_COLLAPSE_PAIR) { + std::swap(edge, pair); + } + p_collapsing_verts(edge, pair, &newv, &keepv); + + if (!(newv->flag & PVERT_PIN)) { + newv->uv[0] = keepv->uv[0]; + newv->uv[1] = keepv->uv[1]; + } + } +} + +static bool p_validate_corrected_coords_point(const PEdge *corr_e, + const float corr_co1[3], + const float corr_co2[3], + const float min_area, + const float min_angle_cos) +{ + /* Check whether the given corrected coordinates don't result in any other triangle with area + * lower than min_area. */ + + const PVert *corr_v = corr_e->vert; + const PEdge *e = corr_v->edge; + + do { + float area; + + if (e == corr_e) { + continue; + } + + if (!(e->face->flag & PFACE_DONE) && (e != corr_e)) { + continue; + } + + if (e->next->next == corr_e->pair) { + PVert *other_v = e->next->next->vert; + area = area_tri_v3(corr_co1, corr_co2, other_v->co); + } + else { + const PVert *other_v1 = e->next->vert; + const PVert *other_v2 = e->next->next->vert; + area = area_tri_v3(corr_co1, other_v1->co, other_v2->co); + } + + if (area < min_area) { + return false; + } + + float f_cos[3]; + + if (e->next->next == corr_e->pair) { + PVert *other_v = e->next->next->vert; + p_triangle_cos(corr_co1, corr_co2, other_v->co, f_cos, f_cos + 1, f_cos + 2); + } + else { + const PVert *other_v1 = e->next->vert; + const PVert *other_v2 = e->next->next->vert; + p_triangle_cos(corr_co1, other_v1->co, other_v2->co, f_cos, f_cos + 1, f_cos + 2); + } + + for (int i = 0; i < 3; i++) { + if (f_cos[i] > min_angle_cos) { + return false; + } + } + + } while ((e = p_wheel_edge_next(e)) && (e != corr_v->edge)); + + return true; +} + +static bool p_validate_corrected_coords(const PEdge *corr_e, + const float corr_co[3], + float min_area, + float min_angle_cos) +{ + /* Check whether the given corrected coordinates don't result in any other triangle with area + * lower than min_area. */ + + const PVert *corr_v = corr_e->vert; + const PEdge *e = corr_v->edge; + + do { + if (!(e->face->flag & PFACE_DONE) && (e != corr_e)) { + continue; + } + + const PVert *other_v1 = e->next->vert; + const PVert *other_v2 = e->next->next->vert; + + const float area = area_tri_v3(corr_co, other_v1->co, other_v2->co); + + if (area < min_area) { + return false; + } + + float f_cos[3]; + p_triangle_cos(corr_co, other_v1->co, other_v2->co, f_cos, f_cos + 1, f_cos + 2); + + for (int i = 0; i < 3; i++) { + if (f_cos[i] > min_angle_cos) { + return false; + } + } + + } while ((e = p_wheel_edge_next(e)) && (e != corr_v->edge)); + + return true; +} + +static bool p_edge_matrix(float R[3][3], const float edge_dir[3]) +{ + static const constexpr float eps = 1.0e-5; + static const constexpr float n1[3] = {0.0f, 0.0f, 1.0f}; + static const constexpr float n2[3] = {0.0f, 1.0f, 0.0f}; + + float edge_len = len_v3(edge_dir); + if (edge_len < eps) { + return false; + } + + float edge_dir_norm[3]; + copy_v3_v3(edge_dir_norm, edge_dir); + mul_v3_fl(edge_dir_norm, 1.0f / edge_len); + + float normal_dir[3]; + cross_v3_v3v3(normal_dir, edge_dir_norm, n1); + float normal_len = len_v3(normal_dir); + + if (normal_len < eps) { + cross_v3_v3v3(normal_dir, edge_dir_norm, n2); + normal_len = len_v3(normal_dir); + + if (normal_len < eps) { + return false; + } + } + + mul_v3_fl(normal_dir, 1.0f / normal_len); + + float tangent_dir[3]; + cross_v3_v3v3(tangent_dir, edge_dir_norm, normal_dir); + + R[0][0] = edge_dir_norm[0]; + R[1][0] = edge_dir_norm[1]; + R[2][0] = edge_dir_norm[2]; + + R[0][1] = normal_dir[0]; + R[1][1] = normal_dir[1]; + R[2][1] = normal_dir[2]; + + R[0][2] = tangent_dir[0]; + R[1][2] = tangent_dir[1]; + R[2][2] = tangent_dir[2]; + + return true; +} + +static bool p_edge_matrix(float R[3][3], const PEdge *e) +{ + float edge_dir[3]; + copy_v3_v3(edge_dir, e->next->vert->co); + sub_v3_v3(edge_dir, e->vert->co); + + return p_edge_matrix(R, edge_dir); +} + +static const float CORR_ZERO_AREA_EPS = 1.0e-10f; + +static bool p_chart_correct_degenerate_triangle_point(PFace *f, + float min_area, + float min_angle_cos) +{ + static const float3 ref_edges[] = { + {1.0f, 0.0f, 0.0f}, + {0.0f, 1.0f, 0.0f}, + {0.0f, 0.0f, 1.0f}, + {0.0f, 1.0f, 1.0f}, + {1.0f, 0.0f, 1.0f}, + {1.0f, 1.0f, 0.0f}, + + {0.0f, 0.5f, 1.0f}, + {0.5f, 0.0f, 1.0f}, + {0.5f, 1.0f, 0.0f}, + + {0.0f, 1.0f, 0.5f}, + {1.0f, 0.0f, 0.5f}, + {1.0f, 0.5f, 0.0f}, + }; + static const int ref_edge_count = ARRAY_SIZE(ref_edges); + static const int LEN_MULTIPLIER_COUNT = 3; + bool corr_co_found = false; + + float corr_len = 2.0f * std::sqrt((min_area + CORR_ZERO_AREA_EPS) / 3.0f / std::sqrt(3.0f)); + + for (int m = 0; m < LEN_MULTIPLIER_COUNT; m++) { + for (int re = 0; re < ref_edge_count; re++) { + float M[3][3]; + + if (!p_edge_matrix(M, ref_edges[re])) { + continue; + } + + float corr_co[3][3]; + PEdge *e = f->edge; + + for (int i = 0; i < 3; i++) { + const float angle = (float(i) / 3.0f) * float(2.0 * M_PI); + float corr_dir[3] = {0.0f, cos(angle), sin(angle)}; + + float corr_len_multiplied = corr_len * (m + 1); + + mul_m3_v3(M, corr_dir); + mul_v3_fl(corr_dir, corr_len_multiplied); + + copy_v3_v3(corr_co[i], e->vert->co); + add_v3_v3(corr_co[i], corr_dir); + e = e->next; + } + + e = f->edge; + for (int i = 0; i < 3; i++) { + if (!p_validate_corrected_coords_point( + e, corr_co[i], corr_co[(i + 1) % 3], min_area, min_angle_cos)) + { + return false; + } + + e = e->next; + } + + e = f->edge; + for (int i = 0; i < 3; i++) { + copy_v3_v3(e->vert->co, corr_co[i]); + e = e->next; + } + + corr_co_found = true; + break; + } + + if (corr_co_found) { + break; + } + } + + if (!corr_co_found) { + return false; + } + + f->flag |= PFACE_DONE; + return true; +} + +static bool p_chart_correct_degenerate_triangles2(PChart *chart, float min_area, float min_angle) +{ + static const float eps = 1.0e-6; + + float min_angle_cos = std::cos(min_angle); + float min_angle_sin = std::sin(min_angle + CORR_ZERO_AREA_EPS); + + for (PFace *f = chart->faces; f; f = f->nextlink) { + if (f->flag & PFACE_DONE) { + continue; + } + + float face_area = p_face_area(f); + + PEdge *max_edge = nullptr; + float max_edge_len = -std::numeric_limits::infinity(); + + PEdge *min_edge = nullptr; + float min_edge_len = std::numeric_limits::infinity(); + + PEdge *middle_edge = nullptr; + + PEdge *e = f->edge; + do { + float len = p_edge_length(e); + + if (len > max_edge_len) { + max_edge = e; + max_edge_len = len; + + middle_edge = max_edge->next == min_edge ? min_edge->next : max_edge->next; + } + + if (len < min_edge_len) { + min_edge = e; + min_edge_len = len; + + middle_edge = min_edge->next == max_edge ? max_edge->next : min_edge->next; + } + + e = e->next; + } while (e != f->edge); + + BLI_assert(max_edge); + BLI_assert(min_edge); + BLI_assert(middle_edge); + + bool small_uniside_tri = (face_area <= min_area) && (min_edge == max_edge); + + if ((max_edge_len < eps) || small_uniside_tri) { + p_chart_correct_degenerate_triangle_point(f, min_area, min_angle_cos); + continue; + } + + if (min_edge == max_edge) { + BLI_assert(face_area > min_area); + f->flag |= PFACE_DONE; + continue; + } + + BLI_assert(middle_edge != max_edge); + BLI_assert(middle_edge != min_edge); + + float M[3][3]; + if (!p_edge_matrix(M, max_edge)) { + continue; + } + + float max_face_cos = + middle_edge->next == max_edge ? + p_vec_cos(middle_edge->vert->co, max_edge->vert->co, min_edge->vert->co) : + p_vec_cos(max_edge->vert->co, middle_edge->vert->co, min_edge->vert->co); + + float angle_corr_len = max_face_cos > min_angle_cos ? + p_edge_length(middle_edge) * min_angle_sin : + 0.0f; + + if ((face_area > min_area) && (angle_corr_len == 0.0f)) { + f->flag |= PFACE_DONE; + continue; + } + + float corr_len = (min_area + CORR_ZERO_AREA_EPS) * 2.0f / max_edge_len; + corr_len = std::max(corr_len, angle_corr_len); + + PEdge *corr_e = max_edge->next->next; + PVert *corr_v = corr_e->vert; + + /* Check 4 distinct directions. */ + static const constexpr int DIR_COUNT = 16; + static const constexpr int LEN_MULTIPLIER_COUNT = 2; + float corr_co[3]; + bool corr_co_found = false; + + for (int m = 0; m < LEN_MULTIPLIER_COUNT; m++) { + for (int d = 0; d < DIR_COUNT; d++) { + const float angle = (float(d) / DIR_COUNT) * float(2.0 * M_PI); + float corr_dir[3] = {0.0f, cos(angle), sin(angle)}; + + const float corr_len_multiplied = corr_len * (m + 1); + + mul_m3_v3(M, corr_dir); + mul_v3_fl(corr_dir, corr_len_multiplied); + + copy_v3_v3(corr_co, corr_v->co); + add_v3_v3(corr_co, corr_dir); + + if (p_validate_corrected_coords(corr_e, corr_co, min_area, min_angle_cos)) { + corr_co_found = true; + break; + } + } + + if (corr_co_found) { + break; + } + } + + if (!corr_co_found) { + continue; + } + + f->flag |= PFACE_DONE; + copy_v3_v3(corr_v->co, corr_co); + } + + return true; +} + +#ifndef NDEBUG + +static bool p_validate_triangle_angles(const PVert *vert1, + const PVert *vert2, + const PVert *vert3, + const float min_angle_cos) +{ + float t_cos[3]; + p_triangle_cos(vert1->co, vert2->co, vert3->co, t_cos, t_cos + 1, t_cos + 2); + + for (int i = 0; i < 3; i++) { + if (t_cos[i] > min_angle_cos) { + return false; + } + } + + return true; +} + +#endif + +static bool UNUSED_FUNCTION_NO_SLIM(p_chart_correct_degenerate_triangles)(PChart *chart, + const float min_area, + const float min_angle) +{ + /* Look for degenerate triangles: triangles with angles lower than `min_angle` or having area + * lower than `min_area` and try to correct vertex coordinates so that the resulting triangle is + * not degenerate. + * + * The return value indicates whether all triangles could be corrected. + */ + + bool ret = p_chart_correct_degenerate_triangles2(chart, min_area, min_angle); + +#ifndef NDEBUG + float min_angle_cos = std::cos(min_angle - CORR_ZERO_AREA_EPS); +#endif + + for (PFace *f = chart->faces; f; f = f->nextlink) { + if (!(f->flag & PFACE_DONE)) { + ret = false; + } + else { +#ifndef NDEBUG + float f_area = p_face_area(f); + BLI_assert(f_area > (min_area - CORR_ZERO_AREA_EPS)); + + PVert *vert1 = f->edge->vert; + PVert *vert2 = f->edge->next->vert; + PVert *vert3 = f->edge->next->next->vert; + + BLI_assert(p_validate_triangle_angles(vert1, vert2, vert3, min_angle_cos)); +#endif + } + + f->flag &= ~PFACE_DONE; + } + + return ret; +} + +/** \} */ + +/* -------------------------------------------------------------------- */ +/** \name SLIM Implementation (Private API) + * \{ */ + +#ifdef WITH_UV_SLIM + +/** + * Get SLIM parameters from the scene. + */ +static slim::MatrixTransfer *slim_matrix_transfer(const ParamSlimOptions *slim_options) +{ + slim::MatrixTransfer *mt = new slim::MatrixTransfer(); + + mt->use_weights = slim_options->weight_influence > 0.0f; + mt->weight_influence = slim_options->weight_influence; + mt->n_iterations = slim_options->iterations; + mt->reflection_mode = slim_options->no_flip; + mt->skip_initialization = slim_options->skip_init; + + return mt; +} + +/** + * For one chart, allocate memory. If no accurate estimate + * (e.g. for number of pinned vertices) overestimate and correct later. + */ +static void slim_allocate_matrices(const PChart *chart, slim::MatrixTransferChart *mt_chart) +{ + mt_chart->verts_num = chart->nverts; + mt_chart->faces_num = chart->nfaces; + mt_chart->edges_num = chart->nedges; + + mt_chart->v_matrices.resize(mt_chart->verts_num * 3); + mt_chart->uv_matrices.resize(mt_chart->verts_num * 2); + mt_chart->pp_matrices.resize(mt_chart->verts_num * 2); + + mt_chart->f_matrices.resize(mt_chart->faces_num * 3); + mt_chart->p_matrices.resize(mt_chart->verts_num); + mt_chart->b_vectors.resize(mt_chart->verts_num); + /* Also clear memory for weight vectors, hence 'new' followed by `()`. */ + mt_chart->w_vectors.resize(mt_chart->verts_num, 0.0); + + mt_chart->e_matrices.resize(mt_chart->edges_num * 2 * 2); + mt_chart->el_vectors.resize(mt_chart->edges_num * 2); +} + +/** + * Transfer edges and edge lengths. + */ +static void slim_transfer_edges(PChart *chart, slim::MatrixTransferChart *mt_chart) +{ + std::vector &E = mt_chart->e_matrices; + std::vector &EL = mt_chart->el_vectors; + + PEdge *outer; + p_chart_boundaries(chart, &outer); + + PEdge *be = outer; + int eid = 0; + + static const float DOUBLED_VERT_THRESHOLD = 1.0e-5; + + do { + E[eid] = be->vert->slim_id; + const float edge_len = p_edge_length(be); + EL[eid] = edge_len; + + /* Temporary solution : SLIM doesn't support doubled vertices for now. */ + if (edge_len < DOUBLED_VERT_THRESHOLD) { + mt_chart->succeeded = false; + } + + be = p_boundary_edge_next(be); + E[eid + mt_chart->edges_num + mt_chart->boundary_vertices_num] = be->vert->slim_id; + eid++; + + } while (be != outer); + + for (PEdge *e = chart->edges; e; e = e->nextlink) { + PEdge *e1 = e->next; + + E[eid] = e->vert->slim_id; + const float edge_len = p_edge_length(e); + EL[eid] = edge_len; + + /* Temporary solution : SLIM doesn't support doubled vertices for now. */ + if (edge_len < DOUBLED_VERT_THRESHOLD) { + mt_chart->succeeded = false; + } + + E[eid + mt_chart->edges_num + mt_chart->boundary_vertices_num] = e1->vert->slim_id; + eid++; + } +} + +/** + * Transfer vertices and pinned information. + */ +static void slim_transfer_vertices(const PChart *chart, + slim::MatrixTransferChart *mt_chart, + slim::MatrixTransfer *mt) +{ + int r = mt_chart->verts_num; + std::vector &v_mat = mt_chart->v_matrices; + std::vector &uv_mat = mt_chart->uv_matrices; + std::vector &p_mat = mt_chart->p_matrices; + std::vector &pp_mat = mt_chart->pp_matrices; + std::vector &w_vec = mt_chart->w_vectors; + + int p_vid = 0; + int vid = mt_chart->boundary_vertices_num; + + /* For every vertex, fill up V matrix and P matrix (pinned vertices) */ + for (PVert *v = chart->verts; v; v = v->nextlink) { + if (!v->on_boundary_flag) { + if (mt->use_weights) { + w_vec[vid] = v->weight; + } + + v->slim_id = vid; + vid++; + } + + v_mat[v->slim_id] = v->co[0]; + v_mat[r + v->slim_id] = v->co[1]; + v_mat[2 * r + v->slim_id] = v->co[2]; + + uv_mat[v->slim_id] = v->uv[0]; + uv_mat[r + v->slim_id] = v->uv[1]; + + if (v->flag & PVERT_PIN || (mt->is_minimize_stretch && !(v->flag & PVERT_SELECT))) { + mt_chart->pinned_vertices_num += 1; + + p_mat[p_vid] = v->slim_id; + pp_mat[2 * p_vid] = double(v->uv[0]); + pp_mat[2 * p_vid + 1] = double(v->uv[1]); + p_vid += 1; + } + } +} + +/** + * Transfer boundary vertices. + */ +static void slim_transfer_boundary_vertices(PChart *chart, + slim::MatrixTransferChart *mt_chart, + const slim::MatrixTransfer *mt) +{ + std::vector &b_vec = mt_chart->b_vectors; + std::vector &w_vec = mt_chart->w_vectors; + + /* For every vertex, set slim_flag to 0 */ + for (PVert *v = chart->verts; v; v = v->nextlink) { + v->on_boundary_flag = false; + } + + /* Find vertices on boundary and save into vector B */ + PEdge *outer; + p_chart_boundaries(chart, &outer); + + PEdge *be = outer; + int vid = 0; + + do { + if (mt->use_weights) { + w_vec[vid] = be->vert->weight; + } + + mt_chart->boundary_vertices_num += 1; + be->vert->slim_id = vid; + be->vert->on_boundary_flag = true; + b_vec[vid] = vid; + + vid += 1; + be = p_boundary_edge_next(be); + + } while (be != outer); +} + +/** + * Transfer faces. + */ +static void slim_transfer_faces(const PChart *chart, slim::MatrixTransferChart *mt_chart) +{ + int fid = 0; + + for (PFace *f = chart->faces; f; f = f->nextlink) { + PEdge *e = f->edge; + PEdge *e1 = e->next; + PEdge *e2 = e1->next; + + int r = mt_chart->faces_num; + std::vector &F = mt_chart->f_matrices; + + F[fid] = e->vert->slim_id; + F[r + fid] = e1->vert->slim_id; + F[2 * r + fid] = e2->vert->slim_id; + fid++; + } +} + +/** + * Conversion Function to build matrix for SLIM Parametrization. + */ +static void slim_convert_blender(ParamHandle *phandle, slim::MatrixTransfer *mt) +{ + static const float SLIM_CORR_MIN_AREA = 1.0e-8; + static const float SLIM_CORR_MIN_ANGLE = DEG2RADF(1.0); + + mt->charts.resize(phandle->ncharts); + + for (int i = 0; i < phandle->ncharts; i++) { + PChart *chart = phandle->charts[i]; + slim::MatrixTransferChart *mt_chart = &mt->charts[i]; + + p_chart_correct_degenerate_triangles(chart, SLIM_CORR_MIN_AREA, SLIM_CORR_MIN_ANGLE); + + mt_chart->succeeded = true; + mt_chart->pinned_vertices_num = 0; + mt_chart->boundary_vertices_num = 0; + + /* Allocate memory for matrices of Vertices,Faces etc. for each chart. */ + slim_allocate_matrices(chart, mt_chart); + + /* For each chart, fill up matrices. */ + slim_transfer_boundary_vertices(chart, mt_chart, mt); + slim_transfer_vertices(chart, mt_chart, mt); + slim_transfer_edges(chart, mt_chart); + slim_transfer_faces(chart, mt_chart); + + mt_chart->pp_matrices.resize(mt_chart->pinned_vertices_num * 2); + mt_chart->pp_matrices.shrink_to_fit(); + + mt_chart->p_matrices.resize(mt_chart->pinned_vertices_num); + mt_chart->p_matrices.shrink_to_fit(); + + mt_chart->b_vectors.resize(mt_chart->boundary_vertices_num); + mt_chart->b_vectors.shrink_to_fit(); + + mt_chart->e_matrices.resize((mt_chart->edges_num + mt_chart->boundary_vertices_num) * 2); + mt_chart->e_matrices.shrink_to_fit(); + } +} + +static void slim_transfer_data_to_slim(ParamHandle *phandle, const ParamSlimOptions *slim_options) +{ + slim::MatrixTransfer *mt = slim_matrix_transfer(slim_options); + + slim_convert_blender(phandle, mt); + phandle->slim_mt = mt; +} + +/** + * Set UV on each vertex after SLIM parametrization, for each chart. + */ +static void slim_flush_uvs(ParamHandle *phandle, + slim::MatrixTransfer *mt, + int *count_changed, + int *count_failed, + bool live = false) +{ + int vid; + PVert *v; + + for (int i = 0; i < phandle->ncharts; i++) { + PChart *chart = phandle->charts[i]; + slim::MatrixTransferChart *mt_chart = &mt->charts[i]; + + if (mt_chart->succeeded) { + if (count_changed) { + (*count_changed)++; + } + + std::vector &UV = mt_chart->uv_matrices; + for (v = chart->verts; v; v = v->nextlink) { + if (v->flag & PVERT_PIN) { + continue; + } + + vid = v->slim_id; + v->uv[0] = UV[vid]; + v->uv[1] = UV[mt_chart->verts_num + vid]; + } + } + else { + if (count_failed) { + (*count_failed)++; + } + + if (!live) { + for (v = chart->verts; v; v = v->nextlink) { + v->uv[0] = 0.0f; + v->uv[1] = 0.0f; + } + } + } + } +} + +/** + * Cleanup memory. + */ +static void slim_free_matrix_transfer(ParamHandle *phandle) +{ + delete phandle->slim_mt; + phandle->slim_mt = nullptr; +} + +static void slim_get_pinned_vertex_data(ParamHandle *phandle, + PChart *chart, + slim::MatrixTransferChart &mt_chart, + slim::PinnedVertexData &pinned_vertex_data) +{ + std::vector &pinned_vertex_indices = pinned_vertex_data.pinned_vertex_indices; + std::vector &pinned_vertex_positions_2D = pinned_vertex_data.pinned_vertex_positions_2D; + std::vector &selected_pins = pinned_vertex_data.selected_pins; + + pinned_vertex_indices.clear(); + pinned_vertex_positions_2D.clear(); + selected_pins.clear(); + + /* Boundary vertices have lower slim_ids, process them first. */ + PEdge *outer; + p_chart_boundaries(chart, &outer); + PEdge *be = outer; + do { + if (be->vert->flag & PVERT_PIN) { + /* Reload vertex position. */ + p_vert_load_pin_select_uvs(phandle, be->vert); + + if (be->vert->flag & PVERT_SELECT) { + selected_pins.push_back(be->vert->slim_id); + } + + pinned_vertex_indices.push_back(be->vert->slim_id); + pinned_vertex_positions_2D.push_back(be->vert->uv[0]); + pinned_vertex_positions_2D.push_back(be->vert->uv[1]); + } + be = p_boundary_edge_next(be); + } while (be != outer); + + PVert *v; + for (v = chart->verts; v; v = v->nextlink) { + if (!v->on_boundary_flag && (v->flag & PVERT_PIN)) { + /* Reload `v`. */ + p_vert_load_pin_select_uvs(phandle, v); + + if (v->flag & PVERT_SELECT) { + selected_pins.push_back(v->slim_id); + } + pinned_vertex_indices.push_back(v->slim_id); + pinned_vertex_positions_2D.push_back(v->uv[0]); + pinned_vertex_positions_2D.push_back(v->uv[1]); + } + } + + mt_chart.pinned_vertices_num = pinned_vertex_indices.size(); +} + +#endif /* WITH_UV_SLIM */ + +/** \} */ + +/* -------------------------------------------------------------------- */ +/** \name SLIM Integration (Public API) + * \{ */ + +void uv_parametrizer_slim_solve(ParamHandle *phandle, + const ParamSlimOptions *slim_options, + int *count_changed, + int *count_failed) +{ +#ifdef WITH_UV_SLIM + slim_transfer_data_to_slim(phandle, slim_options); + slim::MatrixTransfer *mt = phandle->slim_mt; + + mt->parametrize(); + + slim_flush_uvs(phandle, mt, count_changed, count_failed); + slim_free_matrix_transfer(phandle); +#else + *count_changed = 0; + *count_failed = 0; + UNUSED_VARS(phandle, slim_options, count_changed, count_failed); +#endif /* !WITH_UV_SLIM */ +} + +void uv_parametrizer_slim_live_begin(ParamHandle *phandle, const ParamSlimOptions *slim_options) +{ +#ifdef WITH_UV_SLIM + slim_transfer_data_to_slim(phandle, slim_options); + slim::MatrixTransfer *mt = phandle->slim_mt; + + for (int i = 0; i < phandle->ncharts; i++) { + for (PFace *f = phandle->charts[i]->faces; f; f = f->nextlink) { + p_face_backup_uvs(f); + } + } + + for (int i = 0; i < phandle->ncharts; i++) { + PChart *chart = phandle->charts[i]; + slim::MatrixTransferChart &mt_chart = mt->charts[i]; + + bool select = false, deselect = false; + + /* Give vertices matrix indices and count pins. */ + for (PVert *v = chart->verts; v; v = v->nextlink) { + if (v->flag & PVERT_PIN) { + if (v->flag & PVERT_SELECT) { + select = true; + } + } + + if (!(v->flag & PVERT_SELECT)) { + deselect = true; + } + } + + if (!select || !deselect) { + chart->skip_flush = true; + mt_chart.succeeded = false; + continue; + } + + mt->setup_slim_data(mt_chart); + } +#else + UNUSED_VARS(phandle, slim_options); +#endif /* !WITH_UV_SLIM */ +} + +void uv_parametrizer_slim_stretch_iteration(ParamHandle *phandle, const float blend) +{ +#ifdef WITH_UV_SLIM + slim::MatrixTransfer *mt = phandle->slim_mt; + + /* Do one iteration and transfer UVs. */ + for (int i = 0; i < phandle->ncharts; i++) { + mt->charts[i].parametrize_single_iteration(); + mt->charts[i].transfer_uvs_blended(blend); + } + + /* Assign new UVs back to each vertex. */ + slim_flush_uvs(phandle, mt, nullptr, nullptr); +#else + UNUSED_VARS(phandle, blend); +#endif /* !WITH_UV_SLIM */ +} + +void uv_parametrizer_slim_live_solve_iteration(ParamHandle *phandle) +{ +#ifdef WITH_UV_SLIM + slim::MatrixTransfer *mt = phandle->slim_mt; + + /* Do one iteration and transfer UVs */ + for (int i = 0; i < phandle->ncharts; i++) { + PChart *chart = phandle->charts[i]; + slim::MatrixTransferChart &mt_chart = mt->charts[i]; + + if (!mt_chart.data) { + continue; + } + + slim_get_pinned_vertex_data(phandle, chart, mt_chart, mt->pinned_vertex_data); + + mt->parametrize_live(mt_chart, mt->pinned_vertex_data); + } + + /* Assign new UVs back to each vertex. */ + const bool live = true; + slim_flush_uvs(phandle, mt, nullptr, nullptr, live); +#else + UNUSED_VARS(phandle); +#endif /* !WITH_UV_SLIM */ +} + +void uv_parametrizer_slim_live_end(ParamHandle *phandle) +{ +#ifdef WITH_UV_SLIM + slim::MatrixTransfer *mt = phandle->slim_mt; + + for (int i = 0; i < phandle->ncharts; i++) { + slim::MatrixTransferChart *mt_chart = &mt->charts[i]; + mt_chart->free_slim_data(); + } + + slim_free_matrix_transfer(phandle); +#else + UNUSED_VARS(phandle); +#endif /* WITH_UV_SLIM */ +} + +bool uv_parametrizer_is_slim(const ParamHandle *phandle) +{ +#ifdef WITH_UV_SLIM + return phandle->slim_mt != nullptr; +#else + UNUSED_VARS(phandle); + return false; +#endif +} + +/** \} */ + } // namespace blender::geometry diff --git a/source/blender/makesdna/DNA_scene_defaults.h b/source/blender/makesdna/DNA_scene_defaults.h index 67bb41751f2..9087bc6b9f1 100644 --- a/source/blender/makesdna/DNA_scene_defaults.h +++ b/source/blender/makesdna/DNA_scene_defaults.h @@ -368,9 +368,15 @@ .object_flag = SCE_OBJECT_MODE_LOCK, \ .doublimit = 0.001, \ .vgroup_weight = 1.0f, \ + \ .uvcalc_margin = 0.001f, \ .uvcalc_flag = UVCALC_TRANSFORM_CORRECT_SLIDE, \ .unwrapper = UVCALC_UNWRAP_METHOD_CONFORMAL, \ + .uvcalc_iterations = 10, \ + /* See struct member doc-string regarding this name. */ \ + .uvcalc_weight_group = "uv_importance", \ + .uvcalc_weight_factor = 1.0, \ + \ .select_thresh = 0.01f, \ \ .selectmode = SCE_SELECT_VERTEX, \ diff --git a/source/blender/makesdna/DNA_scene_types.h b/source/blender/makesdna/DNA_scene_types.h index 5bd5742b279..93cf9075923 100644 --- a/source/blender/makesdna/DNA_scene_types.h +++ b/source/blender/makesdna/DNA_scene_types.h @@ -1602,6 +1602,20 @@ typedef struct ToolSettings { float uvcalc_margin; + int uvcalc_iterations; + float uvcalc_weight_factor; + + /** + * Regarding having a single vertex group for all meshes. + * In most cases there is no expectation for the names used for vertex groups. + * UV weights is a fairly specific feature for unwrapping and in this case + * users are expected to use the name `uv_importance`. + * While we could support setting a different group per mesh (similar to the active group). + * This isn't all that useful in practice, so use a "default" name instead. + * This approach may be reworked after gathering feedback from users. + */ + char uvcalc_weight_group[64]; /* MAX_VGROUP_NAME */ + /* Auto-IK. */ /** Runtime only. */ short autoik_chainlen; @@ -2695,6 +2709,7 @@ enum { enum { UVCALC_UNWRAP_METHOD_ANGLE = 0, UVCALC_UNWRAP_METHOD_CONFORMAL = 1, + UVCALC_UNWRAP_METHOD_MINIMUM_STRETCH = 2, }; /** #ToolSettings::uvcalc_flag */ @@ -2710,6 +2725,10 @@ enum { UVCALC_TRANSFORM_CORRECT = 1 << 4, /** Keep equal values merged while correcting custom-data. */ UVCALC_TRANSFORM_CORRECT_KEEP_CONNECTED = 1 << 5, + /** Prevent unwrap that flips. */ + UVCALC_UNWRAP_NO_FLIP = 1 << 6, + /** Use importance weights. */ + UVCALC_UNWRAP_USE_WEIGHTS = 1 << 7, }; /** #ToolSettings::uv_flag */ diff --git a/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc b/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc index c6668225e65..cfb7b0f32a1 100644 --- a/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc +++ b/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc @@ -73,6 +73,7 @@ static VArray construct_uv_gvarray(const Mesh &mesh, mp_vkeys.data(), mp_co.data(), mp_uv.data(), + nullptr, mp_pin.data(), mp_select.data()); }); diff --git a/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc b/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc index fcf5191c1ac..b2eef4f15ab 100644 --- a/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc +++ b/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc @@ -104,6 +104,7 @@ static VArray construct_uv_gvarray(const Mesh &mesh, mp_vkeys.data(), mp_co.data(), mp_uv.data(), + nullptr, mp_pin.data(), mp_select.data()); });