UV: add support for the SLIM unwrapping algorithm

Integrate an existing implementation of the SLIM unwrapping algorithm
into Blender. More info about SLIM here:
https://igl.ethz.ch/projects/slim/

This commit is based on the integration code written by Aurel Gruber
for Blender 2.7x (unfinished and never merged with the main branch).

This commit is based on Aurel's code, rebased and further improved.

Details:

- Unwrap has been moved into a sub-menu,
  slim unwrapping is exposed as: "Minimum Stretch".
- Live unwrap with SLIM refines the solutions using a timer.
- When using SLIM there are options to:
  - Set the number of iterations.
  - Weight the influence using vertex weights.
- SLIM can be disabled using the `WITH_UV_SLIM` build option.

Co-authored-by: Aurel Gruber <aurel.gruber@infix.ch>

Ref !114545

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)

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)

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)

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)

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 */
 

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.

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()

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)

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 <Eigen/Dense>
+
+#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<typename VertexPositionType, typename FaceIndicesType>
+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

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 <Eigen/Dense>
+
+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

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 <vector>
+
+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<typename DerivedV, typename DerivedF, typename DerivedC>
+static inline void cotmatrix_entries(const Eigen::PlainObjectBase<DerivedV> &V,
+                                     const Eigen::PlainObjectBase<DerivedF> &F,
+                                     Eigen::PlainObjectBase<DerivedC> &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<typename DerivedC::Scalar, Dynamic, 3> l2;
+  squared_edge_lengths(V, F, l2);
+  /* Compute Edge lenghts. */
+  Matrix<typename DerivedC::Scalar, Dynamic, 3> l;
+  l = l2.array().sqrt();
+
+  /* Double area. */
+  Matrix<typename DerivedC::Scalar, Dynamic, 1> 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<typename DerivedV, typename DerivedF, typename Scalar>
+inline void cotmatrix(const Eigen::PlainObjectBase<DerivedV> &V,
+                      const Eigen::PlainObjectBase<DerivedF> &F,
+                      Eigen::SparseMatrix<Scalar> &L)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  L.resize(V.rows(), V.rows());
+  Matrix<int, Dynamic, 2> 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<Scalar, Dynamic, Dynamic> C;
+  cotmatrix_entries(V, F, C);
+
+  vector<Triplet<Scalar>> 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<Scalar>(source, dest, C(i, e)));
+      IJV.push_back(Triplet<Scalar>(dest, source, C(i, e)));
+      IJV.push_back(Triplet<Scalar>(source, source, -C(i, e)));
+      IJV.push_back(Triplet<Scalar>(dest, dest, -C(i, e)));
+    }
+  }
+  L.setFromTriplets(IJV.begin(), IJV.end());
+}
+
+}  // namespace slim

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 <Eigen/Dense>
+#include <Eigen/Sparse>
+
+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<typename DerivedV, typename DerivedF, typename Scalar>
+inline void cotmatrix(const Eigen::PlainObjectBase<DerivedV> &V,
+                      const Eigen::PlainObjectBase<DerivedF> &F,
+                      Eigen::SparseMatrix<Scalar> &L);
+
+}  // namespace slim
+
+#include "cotmatrix.cpp"

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 <cassert>
+
+#include <BLI_task.hh>
+
+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<typename DerivedX, typename DerivedY, typename DerivedIX>
+static inline void doublearea_sort3(const Eigen::PlainObjectBase<DerivedX> &X,
+                                    const int dim,
+                                    const bool ascending,
+                                    Eigen::PlainObjectBase<DerivedY> &Y,
+                                    Eigen::PlainObjectBase<DerivedIX> &IX)
+{
+  using namespace Eigen;
+  using namespace std;
+  typedef typename Eigen::PlainObjectBase<DerivedY>::Scalar YScalar;
+  Y = X.template cast<YScalar>();
+  /* 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<DerivedIX>::Scalar Index;
+  IX.resize(X.rows(), X.cols());
+  if (dim == 1) {
+    IX.row(0).setConstant(0); /* = Eigen::PlainObjectBase<DerivedIX>::Zero(1,IX.cols());. */
+    IX.row(1).setConstant(1); /* = Eigen::PlainObjectBase<DerivedIX>::Ones (1,IX.cols());. */
+    IX.row(2).setConstant(2); /* = Eigen::PlainObjectBase<DerivedIX>::Ones (1,IX.cols());. */
+  }
+  else {
+    IX.col(0).setConstant(0); /* = Eigen::PlainObjectBase<DerivedIX>::Zero(IX.rows(),1);. */
+    IX.col(1).setConstant(1); /* = Eigen::PlainObjectBase<DerivedIX>::Ones (IX.rows(),1);. */
+    IX.col(2).setConstant(2); /* = Eigen::PlainObjectBase<DerivedIX>::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<typename DerivedV, typename DerivedF, typename DeriveddblA>
+inline void doublearea(const Eigen::PlainObjectBase<DerivedV> &V,
+                       const Eigen::PlainObjectBase<DerivedF> &F,
+                       Eigen::PlainObjectBase<DeriveddblA> &dblA)
+{
+  const int dim = V.cols();
+  /* Only support triangles. */
+  assert(F.cols() == 3);
+  const size_t m = F.rows();
+  /* Compute edge lengths. */
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> 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<DeriveddblA>::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<typename Derivedl, typename DeriveddblA>
+inline void doublearea(const Eigen::PlainObjectBase<Derivedl> &ul,
+                       Eigen::PlainObjectBase<DeriveddblA> &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<typename Derivedl::Scalar, Eigen::Dynamic, 3> 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

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 <Eigen/Dense>
+
+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<typename DerivedV, typename DerivedF, typename DeriveddblA>
+inline void doublearea(const Eigen::PlainObjectBase<DerivedV> &V,
+                       const Eigen::PlainObjectBase<DerivedF> &F,
+                       Eigen::PlainObjectBase<DeriveddblA> &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<typename Derivedl, typename DeriveddblA>
+inline void doublearea(const Eigen::PlainObjectBase<Derivedl> &l,
+                       Eigen::PlainObjectBase<DeriveddblA> &dblA);
+
+}  // namespace slim
+
+#include "doublearea.cpp"

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<typename DerivedV, typename DerivedF, typename DerivedL>
+inline void edge_lengths(const Eigen::PlainObjectBase<DerivedV> &V,
+                         const Eigen::PlainObjectBase<DerivedF> &F,
+                         Eigen::PlainObjectBase<DerivedL> &L)
+{
+  squared_edge_lengths(V, F, L);
+  L = L.array().sqrt().eval();
+}
+
+template<typename DerivedV, typename DerivedF, typename DerivedL>
+inline void squared_edge_lengths(const Eigen::PlainObjectBase<DerivedV> &V,
+                                 const Eigen::PlainObjectBase<DerivedF> &F,
+                                 Eigen::PlainObjectBase<DerivedL> &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

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 <Eigen/Dense>
+
+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<typename DerivedV, typename DerivedF, typename DerivedL>
+inline void edge_lengths(const Eigen::PlainObjectBase<DerivedV> &V,
+                         const Eigen::PlainObjectBase<DerivedF> &F,
+                         Eigen::PlainObjectBase<DerivedL> &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<typename DerivedV, typename DerivedF, typename DerivedL>
+inline void squared_edge_lengths(const Eigen::PlainObjectBase<DerivedV> &V,
+                                 const Eigen::PlainObjectBase<DerivedF> &F,
+                                 Eigen::PlainObjectBase<DerivedL> &L);
+}  // namespace slim
+
+#include "edge_lengths.cpp"

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 <Eigen/Dense>
+#include <vector>
+
+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<double(Eigen::MatrixXd &)> 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<double(Eigen::MatrixXd &)> 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

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 <Eigen/Dense>
+
+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<double(Eigen::MatrixXd &)> energy,
+                                        double cur_energy = -1);
+
+}  // namespace slim
+
+#include "flip_avoiding_line_search.cpp"

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 <Eigen/Dense>
+
+#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<VectorXi>::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

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 <stdio.h>
+
+#include <Eigen/Dense>
+
+#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<MatrixXd> vertex_positions3d = Map<MatrixXd>(NULL, 0, 0);
+  Map<MatrixXd> uv_positions2d = Map<MatrixXd>(NULL, 0, 0);
+  MatrixXd positions_of_pinned_vertices2d;
+  Map<Matrix<double, Dynamic, Dynamic, RowMajor>> positions_of_explicitly_pinned_vertices2d =
+      Map<Matrix<double, Dynamic, Dynamic, RowMajor>>(NULL, 0, 0);
+
+  Map<MatrixXi> faces_by_vertexindices = Map<MatrixXi>(NULL, 0, 0);
+  Map<MatrixXi> edges_by_vertexindices = Map<MatrixXi>(NULL, 0, 0);
+  VectorXi pinned_vertex_indices;
+  Map<VectorXi> explicitly_pinned_vertex_indices = Map<VectorXi>(NULL, 0);
+
+  Map<VectorXd> edge_lengths = Map<VectorXd>(NULL, 0);
+  Map<VectorXi> boundary_vertex_indices = Map<VectorXi>(NULL, 0);
+  Map<VectorXf> weights_per_vertex = Map<VectorXf>(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

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 <Eigen/Dense>
+
+#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 &centroid)
+{
+  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<MatrixXd> 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

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 <stdio.h>
+
+namespace slim {
+
+void transform_initialization_if_necessary(SLIMData &slimData);
+
+}

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 <vector>
+
+#include <Eigen/Geometry>
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/SVD>
+#include <Eigen/SparseCholesky>
+
+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<typename DerivedV, typename DerivedF>
+static inline void grad(const Eigen::PlainObjectBase<DerivedV> &V,
+                        const Eigen::PlainObjectBase<DerivedF> &F,
+                        Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                        bool uniform = false)
+{
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> 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<typename DerivedV::Scalar, 1, 3> v32 = V.row(i3) - V.row(i2);
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v13 = V.row(i1) - V.row(i3);
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v21 = V.row(i2) - V.row(i1);
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> 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<typename DerivedV::Scalar, 1, 3> 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<int> rs;
+  rs.reserve(F.rows() * 4 * 3);
+  std::vector<int> cs;
+  cs.reserve(F.rows() * 4 * 3);
+  std::vector<double> 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<Eigen::Triplet<typename DerivedV::Scalar>> triplets;
+  for (int i = 0; i < (int)vs.size(); ++i) {
+    triplets.push_back(Eigen::Triplet<typename DerivedV::Scalar>(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<typename DerivedA,
+         typename DerivedR,
+         typename DerivedT,
+         typename DerivedU,
+         typename DerivedS,
+         typename DerivedV>
+static inline void polar_svd(const Eigen::PlainObjectBase<DerivedA> &A,
+                             Eigen::PlainObjectBase<DerivedR> &R,
+                             Eigen::PlainObjectBase<DerivedT> &T,
+                             Eigen::PlainObjectBase<DerivedU> &U,
+                             Eigen::PlainObjectBase<DerivedS> &S,
+                             Eigen::PlainObjectBase<DerivedV> &V)
+{
+  using namespace std;
+  Eigen::JacobiSVD<DerivedA> 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<double> &D1,
+                                                   Eigen::SparseMatrix<double> &D2)
+{
+  Eigen::SparseMatrix<double> G;
+  grad(V, F, G);
+  Eigen::SparseMatrix<double> Dx = G.block(0, 0, F.rows(), V.rows());
+  Eigen::SparseMatrix<double> Dy = G.block(F.rows(), 0, F.rows(), V.rows());
+  Eigen::SparseMatrix<double> 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<double>();
+  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<double, 2, 2> Mat2;
+    typedef Eigen::Matrix<double, 2, 1> 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<typename DerivedV, typename DerivedF>
+static inline void local_basis(const Eigen::PlainObjectBase<DerivedV> &V,
+                               const Eigen::PlainObjectBase<DerivedF> &F,
+                               Eigen::PlainObjectBase<DerivedV> &B1,
+                               Eigen::PlainObjectBase<DerivedV> &B2,
+                               Eigen::PlainObjectBase<DerivedV> &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<typename DerivedV::Scalar, 1, 3> v1 =
+        (V.row(F(i, 1)) - V.row(F(i, 0))).normalized();
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> t = V.row(F(i, 2)) - V.row(F(i, 0));
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v3 = v1.cross(t).normalized();
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> 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<double> &A)
+{
+  BLI_assert(s.valid);
+  /* Formula (35) in paper. */
+  std::vector<Eigen::Triplet<double>> 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<double>::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<double>(dx_r, dx_c, weight * val * s.W_11(dx_r)));
+      IJV.push_back(Eigen::Triplet<double>(dx_r, s.v_n + dx_c, weight * val * s.W_12(dx_r)));
+
+      IJV.push_back(Eigen::Triplet<double>(2 * s.f_n + dx_r, dx_c, weight * val * s.W_21(dx_r)));
+      IJV.push_back(
+          Eigen::Triplet<double>(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<double>::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<double>(s.f_n + dy_r, dy_c, weight * val * s.W_11(dy_r)));
+      IJV.push_back(
+          Eigen::Triplet<double>(s.f_n + dy_r, s.v_n + dy_c, weight * val * s.W_12(dy_r)));
+
+      IJV.push_back(Eigen::Triplet<double>(3 * s.f_n + dy_r, dy_c, weight * val * s.W_21(dy_r)));
+      IJV.push_back(
+          Eigen::Triplet<double>(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<double> &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<double> &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<double> &L)
+{
+  BLI_assert(s.valid);
+  /* Formula (35) in paper. */
+  Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
+  buildA(s, A);
+
+  Eigen::SparseMatrix<double> At = A.transpose();
+  At.makeCompressed();
+
+  Eigen::SparseMatrix<double> 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<double> 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<double, 2, 2> 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<double, 2, 2> Mat2;
+    typedef Eigen::Matrix<double, 2, 1> 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<double> L;
+  build_linear_system(s, L);
+
+  /* Solve. */
+  Eigen::VectorXd Uc;
+  SimplicialLDLT<Eigen::SparseMatrix<double>> 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<double(Eigen::MatrixXd &)> 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

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 <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include <stdexcept>
+
+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<double> 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

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 <Eigen/Dense>
+
+#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<SLIMDataPtr::element_type>();
+
+  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

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 <stdlib.h>
+#include <string>
+#include <vector>
+
+#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 <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/SparseCholesky>
+
+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

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 <Eigen/SparseLU>
+
+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<double> &vertex_to_face_indices)
+{
+
+  typedef Eigen::Triplet<double> t;
+  std::vector<t> 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<double> &vertex_to_face_indices)
+{
+
+  typedef Eigen::Triplet<double> t;
+  std::vector<t> 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<double> 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<double> aint(n_unknowns, n_unknowns);
+  Eigen::SparseMatrix<double> abnd(n_unknowns, n_knowns);
+  Eigen::VectorXd z(n_knowns);
+
+  std::vector<Eigen::Triplet<double>> int_triplet_vector;
+  std::vector<Eigen::Triplet<double>> 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<double>(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<double>(rowindex, columnindex, -edge_weight));
+      }
+      else {
+        /* Known point on the border. */
+        columnindex = second_vertex;
+        bnd_triplet_vector.push_back(Eigen::Triplet<double>(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<Eigen::SparseMatrix<double>> 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<int> &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<int> flip_idx;
+  get_flips(f, uv, flip_idx);
+
+  return flip_idx.size();
+}
+
+}  // namespace slim

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 <stdio.h>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+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

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 <memory>
+#include <vector>
+
+namespace slim {
+
+struct SLIMData;
+
+typedef std::unique_ptr<SLIMData> 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<double> v_matrices;
+  /** UV positions of vertices (matrix [verts_num x 2] in a linearized form). */
+  std::vector<double> uv_matrices;
+  /** Positions of pinned vertices (matrix [pinned_vertices_num x 2] in a linearized form). */
+  std::vector<double> pp_matrices;
+  /** Edge lengths. */
+  std::vector<double> el_vectors;
+  /** Weights per vertex. */
+  std::vector<float> w_vectors;
+
+  /** Vertex index triplets making up faces (matrix [faces_num x 3] in a linearized form). */
+  std::vector<int> f_matrices;
+  /** Indices of pinned vertices. */
+  std::vector<int> p_matrices;
+  /** Vertex index tuples making up edges (matrix [edges_num x 2] in a linearized form). */
+  std::vector<int> e_matrices;
+  /** Vertex indices of boundary vertices. */
+  std::vector<int> 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<int> pinned_vertex_indices;
+  std::vector<double> pinned_vertex_positions_2D;
+  std::vector<int> 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<MatrixTransferChart> 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

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()
 

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.menu("IMAGE_MT_uvs_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.separator()
 

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");
+    }
+  }
 }

source/blender/editors/geometry/node_group_operator.cc

@@ -1142,6 +1142,7 @@ static Set<std::string> 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");

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` */
-
+const wmTimer *ED_uvedit_live_unwrap_timer();
-void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit);
+/**
+ * \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);
 

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);
     }
   }
 

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);
 

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:

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}")
 

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);
   }

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<ModifierData *>(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<ModifierData *>(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<ParamKey, BM_DEFAULT_NGON_STACK_SIZE> vkeys(efa->len);
   blender::Array<bool, BM_DEFAULT_NGON_STACK_SIZE> pin(efa->len);
   blender::Array<bool, BM_DEFAULT_NGON_STACK_SIZE> select(efa->len);
   blender::Array<const float *, BM_DEFAULT_NGON_STACK_SIZE> co(efa->len);
   blender::Array<float *, BM_DEFAULT_NGON_STACK_SIZE> uv(efa->len);
+  blender::Array<float, BM_DEFAULT_NGON_STACK_SIZE> 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(
+  blender::geometry::uv_parametrizer_face_add(handle,
-      handle, face_index, i, vkeys.data(), co.data(), uv.data(), pin.data(), select.data());
+                                              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<blender::int2> subsurf_edges = subdiv_mesh->edges();
   const blender::OffsetIndices subsurf_facess = subdiv_mesh->faces();
   const blender::Span<int> subsurf_corner_verts = subdiv_mesh->corner_verts();
+  const blender::Span<MDeformVert> subsurf_deform_verts = subdiv_mesh->deform_verts();
 
   const int *origVertIndices = static_cast<const int *>(
       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<wmWindowManager *>(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<wmWindowManager *>(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(
+  if (options->use_slim) {
-      handle, false, scene->toolsettings->unwrapper == UVCALC_UNWRAP_METHOD_ANGLE);
+    uv_parametrizer_slim_solve(handle, &options->slim, r_count_changed, r_count_failed);
-  blender::geometry::uv_parametrizer_lscm_solve(handle, r_count_changed, r_count_failed);
+  }
-  blender::geometry::uv_parametrizer_lscm_end(handle);
+  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<Object *> 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<Object *> 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,
+  ot->prop = RNA_def_enum(
-               "method",
+      ot->srna,
-               method_items,
+      "method",
-               0,
+      method_items,
-               "Method",
+      tool_settings_default->unwrapper,
-               "Unwrapping method (Angle Based usually gives better results than Conformal, while "
+      "Method",
-               "being somewhat slower)");
+      "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);
 }
 
 /** \} */

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

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:
  *

source/blender/geometry/intern/uv_parametrizer.cc

@@ -6,6 +6,9 @@
  * \ingroup eduv
  */
 
+#include <functional>
+#include <vector>
+
 #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);
+  float weight1, weight2, weight3;
-  e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], e2);
+  if (weight) {
-  e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], e3);
+    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<ParamKey> vkeys_sub(pm);
       Array<const float *> co_sub(pm);
       Array<float *> uv_sub(pm);
+      Array<float> weight_sub(pm);
       Array<bool> pin_sub(pm);
       Array<bool> 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<float>::infinity();
+
+    PEdge *min_edge = nullptr;
+    float min_edge_len = std::numeric_limits<float>::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<int> &E = mt_chart->e_matrices;
+  std::vector<double> &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<double> &v_mat = mt_chart->v_matrices;
+  std::vector<double> &uv_mat = mt_chart->uv_matrices;
+  std::vector<int> &p_mat = mt_chart->p_matrices;
+  std::vector<double> &pp_mat = mt_chart->pp_matrices;
+  std::vector<float> &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<int> &b_vec = mt_chart->b_vectors;
+  std::vector<float> &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<int> &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<double> &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<int> &pinned_vertex_indices = pinned_vertex_data.pinned_vertex_indices;
+  std::vector<double> &pinned_vertex_positions_2D = pinned_vertex_data.pinned_vertex_positions_2D;
+  std::vector<int> &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

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, \

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 */

source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc

@@ -73,6 +73,7 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
                                        mp_vkeys.data(),
                                        mp_co.data(),
                                        mp_uv.data(),
+                                       nullptr,
                                        mp_pin.data(),
                                        mp_select.data());
   });

source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc

@@ -104,6 +104,7 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
                                        mp_vkeys.data(),
                                        mp_co.data(),
                                        mp_uv.data(),
+                                       nullptr,
                                        mp_pin.data(),
                                        mp_select.data());
   });