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

@@ -14,6 +14,7 @@ add_subdirectory(opensubdiv)
 add_subdirectory(mikktspace)
 add_subdirectory(eigen)
 add_subdirectory(sky)
+add_subdirectory(slim)
 
 if(WITH_AUDASPACE)
   add_subdirectory(audaspace)

intern/slim/CMakeLists.txt

@@ -0,0 +1,39 @@
+# SPDX-FileCopyrightText: 2023 Blender Authors
+#
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+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}")

intern/slim/intern/area_compensation.cpp

@@ -0,0 +1,89 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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,17 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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,94 @@
+/* SPDX-FileCopyrightText: 2013 Alec Jacobson
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,40 @@
+/* SPDX-FileCopyrightText: 2014 Alec Jacobson, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,196 @@
+/* SPDX-FileCopyrightText: 2013 Alec Jacobson
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,48 @@
+/* SPDX-FileCopyrightText: 2013 Alec Jacobson, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,43 @@
+/* SPDX-FileCopyrightText: 2013 Alec Jacobson, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,58 @@
+/* SPDX-FileCopyrightText: 2013 Alec Jacobson, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,165 @@
+/* SPDX-FileCopyrightText: 2016 Michael Rabinovich, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,42 @@
+/* SPDX-FileCopyrightText: 2016 Michael Rabinovich, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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,222 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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.n_verts),
+      number_of_faces(chart.n_faces),
+      /* `n_edges` in transferred_data accounts for boundary edges only once. */
+      number_of_edges_twice(chart.n_edges + chart.n_boundary_vertices),
+      number_of_boundary_vertices(chart.n_boundary_vertices),
+      number_of_pinned_vertices(chart.n_pinned_vertices),
+      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,65 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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,235 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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,14 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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,803 @@
+/* SPDX-FileCopyrightText: 2016 Michael Rabinovich, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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
+ *   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,118 @@
+/* SPDX-FileCopyrightText: 2016 Michael Rabinovich, 2023 Blender Authors
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * Derived from libigl, a simple c++ geometry processing library. */
+
+#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). */
+
+  /* weightmap 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,47 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include <Eigen/Dense>
+
+#include "slim.h"
+#include "slim_matrix_transfer.h"
+#include "geometry_data_retrieval.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);
+}
+
+/* Transfers all the matrices from the native part and initialises SLIM. */
+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.n_pinned_vertices = 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,162 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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.n_verts;
+
+  for (int i = 0; i < number_of_vertices; i++) {
+    for (int j = 0; j < 2; j++) {
+      uv_coordinate_array[i * 2 + j] = uv(i, j);
+    }
+  }
+}
+
+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];
+    }
+  }
+}
+
+/* Called from the native part during each iteration of interactive parametrisation.
+ * The blend parameter decides the linear blending between the original UV map and the one
+ * optained from the accumulated SLIM iterations so far. */
+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;
+  }
+}
+
+/* Executes a single iteration of SLIM, must follow a proper setup & initialisation. */
+void MatrixTransferChart::parametrize_single_iteration()
+{
+  int number_of_iterations = 1;
+  try_slim_solve(number_of_iterations);
+}
+
+/* Executes slim iterations during live unwrap. needs to provide new selected-pin positions. */
+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,288 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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 n_verts = uv.rows();
+  int n_edges = e.rows();
+
+  Eigen::SparseMatrix<double> vertex_to_angles(n_verts, n_verts * 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 = n_verts - 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 < n_edges; 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 n_boundary_vertices = vertex_positions.rows();
+  double x, y;
+  double angle = 2 * pi / n_boundary_vertices;
+
+  for (int i = 0; i < n_boundary_vertices; 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,59 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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,95 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#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 n_verts = 0;
+  int n_faces = 0;
+  int n_pinned_vertices = 0;
+  int n_boundary_vertices = 0;
+  int n_edges = 0;
+
+  bool succeeded = false;
+
+  /* Vertex positions. */
+  std::vector<double> v_matrices;
+  /* UV positions of vertices. */
+  std::vector<double> uv_matrices;
+  /* Positions of pinned vertices. */
+  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. */
+  std::vector<int> f_matrices;
+  /* Indices of pinned vertices. */
+  std::vector<int> p_matrices;
+  /* Vertex index tuples making up edges. */
+  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);
+  void parametrize_single_iteration();
+
+  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();
+
+  void parametrize_live(MatrixTransferChart &chart,
+                        const PinnedVertexData &pinned_vertex_data);
+
+  void setup_slim_data(MatrixTransferChart &chart) const;
+};
+
+}  // namespace slim

source/blender/blenkernel/BKE_customdata.hh

@@ -51,6 +51,7 @@ struct BMUVOffsets {
   int select_vert;
   int select_edge;
   int pin;
+  int weight;
 };
 
 /* A data type large enough to hold 1 element from any custom-data layer type. */

source/blender/blenloader/intern/versioning_260.cc

@@ -31,6 +31,7 @@
 #include "DNA_node_types.h"
 #include "DNA_object_fluidsim_types.h"
 #include "DNA_object_types.h"
+#include "DNA_scene_defaults.h"
 #include "DNA_screen_types.h"
 #include "DNA_sequence_types.h"
 #include "DNA_space_types.h"
@@ -2002,8 +2003,8 @@ void blo_do_versions_260(FileData *fd, Library * /*lib*/, Main *bmain)
     /* set a unwrapping margin and ABF by default */
     LISTBASE_FOREACH (Scene *, scene, &bmain->scenes) {
       if (scene->toolsettings->uvcalc_margin == 0.0f) {
-        scene->toolsettings->uvcalc_margin = 0.001f;
-        scene->toolsettings->unwrapper = 0;
+        scene->toolsettings->uvcalc_margin = _DNA_DEFAULT_ToolSettings_UVCalc_Margin;
+        scene->toolsettings->unwrapper = _DNA_DEFAULT_ToolSettings_UVCalc_Unwrapper;
       }
     }
   }

source/blender/blenloader/intern/versioning_300.cc

@@ -45,6 +45,7 @@
 #include "DNA_mesh_types.h"
 #include "DNA_modifier_types.h"
 #include "DNA_movieclip_types.h"
+#include "DNA_scene_defaults.h"
 #include "DNA_screen_types.h"
 #include "DNA_sequence_types.h"
 #include "DNA_space_types.h"
@@ -4556,4 +4557,18 @@ 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 = _DNA_DEFAULT_ToolSettings_UVCalc_Iterations;
+      scene->toolsettings->uvcalc_vertex_group_factor =
+          _DNA_DEFAULT_ToolSettings_UVCalc_VertexGroupFactor;
+      scene->toolsettings->uvcalc_relative_scale = _DNA_DEFAULT_ToolSettings_UVCalc_RelativeScale;
+      scene->toolsettings->uvcalc_allow_flips = _DNA_DEFAULT_ToolSettings_UVCalc_AllowFlips;
+      memset(scene->toolsettings->uvcalc_vertex_group,
+             0,
+             sizeof(scene->toolsettings->uvcalc_vertex_group));
+    }
+  }
 }

source/blender/bmesh/intern/bmesh_query_uv.cc

@@ -24,7 +24,7 @@ BMUVOffsets BM_uv_map_get_offsets_from_layer(const BMesh *bm, const int layer)
   using namespace blender::bke;
   const int layer_index = CustomData_get_layer_index_n(&bm->ldata, CD_PROP_FLOAT2, layer);
   if (layer_index == -1) {
-    return {-1, -1, -1, -1};
+    return {-1, -1, -1, -1, -1};
   }
 
   char const *name = bm->ldata.layers[layer_index].name;
@@ -38,6 +38,7 @@ BMUVOffsets BM_uv_map_get_offsets_from_layer(const BMesh *bm, const int layer)
       &bm->ldata, CD_PROP_BOOL, BKE_uv_map_edge_select_name_get(name, buffer));
   offsets.pin = CustomData_get_offset_named(
       &bm->ldata, CD_PROP_BOOL, BKE_uv_map_pin_name_get(name, buffer));
+  offsets.weight = CustomData_get_offset(&bm->vdata, CD_MDEFORMVERT);
 
   return offsets;
 }
@@ -46,7 +47,7 @@ BMUVOffsets BM_uv_map_get_offsets(const BMesh *bm)
 {
   const int layer = CustomData_get_active_layer(&bm->ldata, CD_PROP_FLOAT2);
   if (layer == -1) {
-    return {-1, -1, -1, -1};
+    return {-1, -1, -1, -1, -1};
   }
   return BM_uv_map_get_offsets_from_layer(bm, layer);
 }

source/blender/editors/include/ED_uvedit.hh

@@ -260,9 +260,9 @@ void ED_uvedit_select_sync_flush(const ToolSettings *ts, BMEditMesh *em, bool se
 
 /* `uvedit_unwrap_ops.cc` */
 
-void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit);
+void ED_uvedit_live_unwrap_begin(bContext *C, Scene *scene, Object *obedit);
 void ED_uvedit_live_unwrap_re_solve();
-void ED_uvedit_live_unwrap_end(short cancel);
+void ED_uvedit_live_unwrap_end(bContext *C, short cancel);
 
 void ED_uvedit_live_unwrap(const Scene *scene, blender::Span<Object *> objects);
 void ED_uvedit_add_simple_uvs(Main *bmain, const Scene *scene, Object *ob);

source/blender/editors/sculpt_paint/sculpt_uv.cc

@@ -588,7 +588,7 @@ static void uv_sculpt_stroke_exit(bContext *C, wmOperator *op)
 {
   SpaceImage *sima = CTX_wm_space_image(C);
   if (sima->flag & SI_LIVE_UNWRAP) {
-    ED_uvedit_live_unwrap_end(false);
+    ED_uvedit_live_unwrap_end(C, false);
   }
   UvSculptData *data = static_cast<UvSculptData *>(op->customdata);
   if (data->timer) {
@@ -905,7 +905,7 @@ 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);
+      ED_uvedit_live_unwrap_begin(C, scene, obedit);
     }
   }
 

source/blender/editors/transform/transform.cc

@@ -1006,6 +1006,9 @@ int transformEvent(TransInfo *t, wmOperator *op, const wmEvent *event)
     t->redraw |= TREDRAW_HARD;
     handled = true;
   }
+  else if (event->type == TIMER) {
+    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,7 @@ static void createTransUVs(bContext *C, TransInfo *t)
     }
 
     if (sima->flag & SI_LIVE_UNWRAP) {
-      ED_uvedit_live_unwrap_begin(t->scene, tc->obedit);
+      ED_uvedit_live_unwrap_begin(C, t->scene, tc->obedit);
     }
 
   finally:

source/blender/editors/transform/transform_generics.cc

@@ -802,7 +802,7 @@ void postTrans(bContext *C, TransInfo *t)
     else {
       SpaceImage *sima = static_cast<SpaceImage *>(t->area->spacedata.first);
       if (sima->flag & SI_LIVE_UNWRAP) {
-        ED_uvedit_live_unwrap_end(t->state == TRANS_CANCEL);
+        ED_uvedit_live_unwrap_end(C, t->state == TRANS_CANCEL);
       }
     }
   }

source/blender/editors/uvedit/CMakeLists.txt

@@ -12,6 +12,7 @@ set(INC
   ../../makesrna
   ../../windowmanager
   ../../../../intern/eigen
+  ../../../../intern/slim
   # RNA_prototypes.h
   ${CMAKE_BINARY_DIR}/source/blender/makesrna
 )
@@ -39,6 +40,7 @@ set(SRC
 set(LIB
   PRIVATE bf::blenlib
   bf_bmesh
+  bf_intern_slim
   PRIVATE bf::depsgraph
   PRIVATE bf::dna
   PRIVATE bf::intern::guardedalloc

source/blender/editors/uvedit/uvedit_ops.cc

@@ -181,9 +181,9 @@ 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(NULL, scene, obedit);
     ED_uvedit_live_unwrap_re_solve();
-    ED_uvedit_live_unwrap_end(0);
+    ED_uvedit_live_unwrap_end(NULL, 0);
   }
 }
 

source/blender/editors/uvedit/uvedit_unwrap_ops.cc

@@ -6,14 +6,17 @@
  * \ingroup eduv
  */
 
+#include <cassert>
 #include <cmath>
 #include <cstdlib>
 #include <cstring>
 
 #include "MEM_guardedalloc.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,6 +39,7 @@
 
 #include "BKE_context.hh"
 #include "BKE_customdata.hh"
+#include "BKE_deform.hh"
 #include "BKE_editmesh.hh"
 #include "BKE_image.h"
 #include "BKE_layer.hh"
@@ -74,31 +78,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 +187,14 @@ 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;
+
+  ParamSlimOptions slim_options;
+  char slim_vertex_group[MAX_VGROUP_NAME];
 };
 
 void blender::geometry::UVPackIsland_Params::setFromUnwrapOptions(const UnwrapOptions &options)
@@ -213,6 +206,153 @@ 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);
+
+  /* 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 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_options.skip_initialization = 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;
+
+    STRNCPY(options.slim_vertex_group, ts->uvcalc_vertex_group);
+
+    options.slim_options.weight_influence = strlen(options.slim_vertex_group) ?
+                                                ts->uvcalc_vertex_group_factor :
+                                                0.0f;
+    options.slim_options.iterations = ts->uvcalc_iterations;
+    options.slim_options.allow_flips = ts->uvcalc_allow_flips;
+  }
+  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");
+
+    RNA_string_get(op->ptr, "vertex_group", options.slim_vertex_group);
+    options.slim_options.weight_influence = strlen(options.slim_vertex_group) ?
+                                                RNA_float_get(op->ptr, "vertex_group_factor") :
+                                                0.0f;
+    options.slim_options.iterations = RNA_int_get(op->ptr, "iterations");
+    options.slim_options.allow_flips = RNA_boolean_get(op->ptr, "allow_flips");
+  }
+
+  options.use_abf = options.method == 0;
+  options.use_slim = options.method == 2;
+
+  /* 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;
+}
+
+static void unwrap_options_sync_flag(
+    wmOperator *op, ToolSettings *ts, const char *property_name, char flag, bool flipped)
+{
+  if (RNA_struct_property_is_set(op->ptr, property_name)) {
+    if (RNA_boolean_get(op->ptr, property_name) ^ flipped) {
+      ts->uvcalc_flag |= flag;
+    }
+    else {
+      ts->uvcalc_flag &= ~flag;
+    }
+  }
+  else {
+    RNA_boolean_set(op->ptr, property_name, ((ts->uvcalc_flag & flag) > 0) ^ flipped);
+  }
+}
+
+static void unwrap_options_sync_toolsettings(wmOperator *op, ToolSettings *ts)
+{
+  /* remember last method for live unwrap */
+  if (RNA_struct_property_is_set(op->ptr, "method")) {
+    ts->unwrapper = RNA_enum_get(op->ptr, "method");
+  }
+  else {
+    RNA_enum_set(op->ptr, "method", ts->unwrapper);
+  }
+
+  /* remember packing margin */
+  if (RNA_struct_property_is_set(op->ptr, "margin")) {
+    ts->uvcalc_margin = RNA_float_get(op->ptr, "margin");
+  }
+  else {
+    RNA_float_set(op->ptr, "margin", ts->uvcalc_margin);
+  }
+
+  if (RNA_struct_property_is_set(op->ptr, "allow_flips")) {
+    ts->uvcalc_allow_flips = RNA_boolean_get(op->ptr, "allow_flips");
+  }
+  else {
+    RNA_boolean_set(op->ptr, "allow_flips", ts->uvcalc_allow_flips);
+  }
+
+  if (RNA_struct_property_is_set(op->ptr, "iterations")) {
+    ts->uvcalc_iterations = RNA_int_get(op->ptr, "iterations");
+  }
+  else {
+    RNA_int_set(op->ptr, "iterations", ts->uvcalc_iterations);
+  }
+
+  if (RNA_struct_property_is_set(op->ptr, "vertex_group_factor")) {
+    ts->uvcalc_vertex_group_factor = RNA_float_get(op->ptr, "vertex_group_factor");
+  }
+  else {
+    RNA_float_set(op->ptr, "vertex_group_factor", ts->uvcalc_vertex_group_factor);
+  }
+
+  if (RNA_struct_property_is_set(op->ptr, "vertex_group")) {
+    RNA_string_get(op->ptr, "vertex_group", ts->uvcalc_vertex_group);
+  }
+  else {
+    RNA_string_set(op->ptr, "vertex_group", ts->uvcalc_vertex_group);
+  }
+
+  unwrap_options_sync_flag(op, ts, "fill_holes", UVCALC_FILLHOLES, false);
+  unwrap_options_sync_flag(op, ts, "correct_aspect", UVCALC_NO_ASPECT_CORRECT, true);
+  unwrap_options_sync_flag(op, ts, "use_subsurf_data", UVCALC_USESUBSURF, false);
+}
+
 static bool uvedit_have_selection(const Scene *scene, BMEditMesh *em, const UnwrapOptions *options)
 {
   BMFace *efa;
@@ -360,13 +500,16 @@ 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_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 +528,26 @@ static void construct_param_handle_face_add(ParamHandle *handle,
     if (options->pin_unselected && !select[i]) {
       pin[i] = true;
     }
+
+    /* optional vertex group weighting */
+    if (offsets.weight >= 0 && cd_weight_index >= 0) {
+      MDeformVert *dv = (MDeformVert *)BM_ELEM_CD_GET_VOID_P(l->v, offsets.weight);
+      weight[i] = BKE_defvert_find_weight(dv, cd_weight_index);
+    }
+    else {
+      weight[i] = 1.0f;
+    }
   }
 
-  blender::geometry::uv_parametrizer_face_add(
-      handle, face_index, i, vkeys.data(), co.data(), uv.data(), pin.data(), select.data());
+  blender::geometry::uv_parametrizer_face_add(handle,
+                                              face_index,
+                                              i,
+                                              vkeys.data(),
+                                              co.data(),
+                                              uv.data(),
+                                              weight.data(),
+                                              pin.data(),
+                                              select.data());
 }
 
 /* Set seams on UV Parametrizer based on options. */
@@ -454,6 +613,8 @@ 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_index = BKE_object_defgroup_name_index(ob, options->slim_vertex_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 +623,7 @@ 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_index);
     }
   }
 
@@ -507,6 +668,8 @@ static ParamHandle *construct_param_handle_multi(const Scene *scene,
       continue;
     }
 
+    const int cd_weight_index = BKE_object_defgroup_name_index(obedit, options->slim_vertex_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 +678,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_index);
       }
     }
 
@@ -608,6 +772,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->slim_vertex_group);
 
   ParamHandle *handle = new blender::geometry::ParamHandle();
 
@@ -630,6 +795,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 +832,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 +863,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 +913,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 +1620,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 +1911,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, nullptr);
   options.topology_from_uvs = true;
   options.only_selected_faces = true;
   options.only_selected_uvs = true;
@@ -1788,36 +1974,43 @@ 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)
+void ED_uvedit_live_unwrap_begin(bContext *C, Scene *scene, Object *obedit)
 {
   ParamHandle *handle = nullptr;
   BMEditMesh *em = BKE_editmesh_from_object(obedit);
-  const bool abf = (scene->toolsettings->unwrapper == 0);
-  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_options.allow_flips = true;
+    options.slim_options.skip_initialization = true;
+    uv_parametrizer_slim_live_begin(handle, &options.slim_options);
+
+    if (C) {
+      BLI_assert(!g_live_unwrap.timer);
+      g_live_unwrap.timer = WM_event_timer_add(CTX_wm_manager(C), CTX_wm_window(C), 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,17 +2032,34 @@ 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]);
     }
   }
 }
 
-void ED_uvedit_live_unwrap_end(short cancel)
+void ED_uvedit_live_unwrap_end(bContext *C, short cancel)
 {
+  if (C && g_live_unwrap.timer) {
+    WM_event_timer_remove(CTX_wm_manager(C), CTX_wm_window(C), g_live_unwrap.timer);
+    g_live_unwrap.timer = NULL;
+  }
+
   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 +2575,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 +2585,14 @@ static void uvedit_unwrap(const Scene *scene,
     handle = construct_param_handle(scene, obedit, em->bm, options, r_count_failed);
   }
 
-  blender::geometry::uv_parametrizer_lscm_begin(
-      handle, false, scene->toolsettings->unwrapper == 0);
-  blender::geometry::uv_parametrizer_lscm_solve(handle, r_count_changed, r_count_failed);
-  blender::geometry::uv_parametrizer_lscm_end(handle);
+  if (options->use_slim) {
+    uv_parametrizer_slim_solve(handle, &options->slim_options, r_count_changed, r_count_failed);
+  }
+  else {
+    blender::geometry::uv_parametrizer_lscm_begin(handle, false, options->use_abf);
+    blender::geometry::uv_parametrizer_lscm_solve(handle, r_count_changed, r_count_failed);
+    blender::geometry::uv_parametrizer_lscm_end(handle);
+  }
 
   blender::geometry::uv_parametrizer_average(handle, true, false, false);
 
@@ -2403,7 +2617,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 +2645,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 +2685,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 +2722,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 +2753,50 @@ 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(NULL, 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") == 2;
+
+  if (is_slim) {
+    uiItemR(col, &ptr, "iterations", UI_ITEM_NONE, nullptr, ICON_NONE);
+    uiItemR(col, &ptr, "allow_flips", UI_ITEM_NONE, nullptr, ICON_NONE);
+
+    col = uiLayoutColumn(layout, true);
+    uiItemR(col, &ptr, "vertex_group", UI_ITEM_NONE, nullptr, ICON_NONE);
+    uiItemR(col, &ptr, "vertex_group_factor", UI_ITEM_NONE, nullptr, ICON_NONE);
+  }
+  else {
+    col = uiLayoutColumn(layout, true);
+    uiItemR(col, &ptr, "fill_holes", UI_ITEM_NONE, nullptr, ICON_NONE);
+  }
+
+  col = uiLayoutColumn(layout, true);
+  uiItemR(col, &ptr, "use_subsurf_data", UI_ITEM_NONE, nullptr, ICON_NONE);
+
+  col = uiLayoutColumn(layout, true);
+  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)
 {
   static const EnumPropertyItem method_items[] = {
       {0, "ANGLE_BASED", 0, "Angle Based", ""},
       {1, "CONFORMAL", 0, "Conformal", ""},
+      {2, "SLIM", 0, "SLIM", ""},
       {0, nullptr, 0, nullptr, nullptr},
   };
 
@@ -2593,23 +2810,26 @@ void UV_OT_unwrap(wmOperatorType *ot)
   ot->exec = unwrap_exec;
   ot->poll = ED_operator_uvmap;
 
+  /* Only draw relevant ui elements */
+  ot->ui = unwrap_draw;
+
   /* properties */
   RNA_def_enum(ot->srna,
                "method",
                method_items,
-               0,
+               _DNA_DEFAULT_ToolSettings_UVCalc_Unwrapper,
                "Method",
                "Unwrapping method (Angle Based usually gives better results than Conformal, while "
                "being somewhat slower)");
   RNA_def_boolean(ot->srna,
                   "fill_holes",
-                  true,
+                  _DNA_DEFAULT_ToolSettings_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,
+                  !(_DNA_DEFAULT_ToolSettings_UVCalc_Flag & UVCALC_NO_ASPECT_CORRECT),
                   "Correct Aspect",
                   "Map UVs taking image aspect ratio into account");
   RNA_def_boolean(
@@ -2626,6 +2846,40 @@ 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,
+                  "allow_flips",
+                  _DNA_DEFAULT_ToolSettings_UVCalc_AllowFlips,
+                  "Allow Flips",
+                  "Allowing flips means that depending on the position of pins, the map may be "
+                  "flipped to lower distortion");
+
+  RNA_def_int(ot->srna,
+              "iterations",
+              _DNA_DEFAULT_ToolSettings_UVCalc_Iterations,
+              0,
+              10000,
+              "Iterations",
+              "Number of Iterations if the SLIM algorithm is used",
+              1,
+              30);
+  RNA_def_string(ot->srna,
+                 "vertex_group",
+                 NULL,
+                 MAX_ID_NAME,
+                 "Vertex Group",
+                 "Vertex group name for modulating the deform");
+  RNA_def_float(
+      ot->srna,
+      "vertex_group_factor",
+      _DNA_DEFAULT_ToolSettings_UVCalc_VertexGroupFactor,
+      -10000.0,
+      10000.0,
+      "Factor",
+      "How much influence the weightmap has for weighted parameterization, 0 being no influence",
+      -10.0,
+      10.0);
 }
 
 /** \} */

source/blender/geometry/CMakeLists.txt

@@ -1,4 +1,4 @@
-# SPDX-FileCopyrightText: 2006 Blender Authors
+# SPDX-FileCopyrightText: 2023 Blender Authors
 #
 # SPDX-License-Identifier: GPL-2.0-or-later
 
@@ -10,6 +10,7 @@ set(INC
   ../functions
   ../makesrna
   ../../../intern/eigen
+  ../../../intern/slim
 )
 
 set(INC_SYS

source/blender/geometry/GEO_uv_parametrizer.hh

@@ -5,6 +5,11 @@
 #pragma once
 
 #include "BLI_sys_types.h" /* for intptr_t support */
+#include "DNA_ID.h"        /* for MAX_ID_NAME */
+
+namespace slim {
+struct MatrixTransfer;
+}
 
 /** \file
  * \ingroup geo
@@ -55,6 +60,9 @@ class ParamHandle {
 
   RNG *rng = nullptr;
   float blend = 0.0f;
+
+  /* SLIM uv unwrapping */
+  slim::MatrixTransfer *slim_mt = nullptr;
 };
 
 /* -------------------------------------------------------------------- */
@@ -84,6 +92,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 +105,35 @@ 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 {
+ public:
+  float weight_influence = 0.0f;
+  int iterations = 0;
+  bool allow_flips = true;
+  bool skip_initialization = false;
+};
+
+void uv_parametrizer_slim_solve(ParamHandle *handle,
+                                const ParamSlimOptions *mt_options,
+                                int *count_changed,
+                                int *count_failed);
+
+void uv_parametrizer_slim_live_begin(ParamHandle *handle, const ParamSlimOptions *mt_options);
+void uv_parametrizer_slim_live_solve_iteration(ParamHandle *handle);
+void uv_parametrizer_slim_live_end(ParamHandle *handle);
+void uv_parametrizer_slim_stretch_iteration(ParamHandle *handle, float blend);
+bool uv_parametrizer_is_slim(ParamHandle *handle);
+
+/** \} */
+
 /* -------------------------------------------------------------------- */
 /** \name Least Squares Conformal Maps:
  *

source/blender/makesdna/DNA_scene_defaults.h

@@ -357,14 +357,29 @@
     .sharp_max = DEG2RADF(180.0f), \
   }
 
+#define _DNA_DEFAULT_ToolSettings_UVCalc_Margin               (0.001f)
+#define _DNA_DEFAULT_ToolSettings_UVCalc_Flag                 (UVCALC_TRANSFORM_CORRECT_SLIDE)
+#define _DNA_DEFAULT_ToolSettings_UVCalc_Unwrapper            (1)
+#define _DNA_DEFAULT_ToolSettings_UVCalc_Iterations           (10)
+#define _DNA_DEFAULT_ToolSettings_UVCalc_VertexGroupFactor    (1.0)
+#define _DNA_DEFAULT_ToolSettings_UVCalc_RelativeScale        (1.0)
+#define _DNA_DEFAULT_ToolSettings_UVCalc_AllowFlips       (1)
+
 #define _DNA_DEFAULT_ToolSettings \
   { \
     .object_flag = SCE_OBJECT_MODE_LOCK, \
     .doublimit = 0.001, \
     .vgroup_weight = 1.0f, \
-    .uvcalc_margin = 0.001f, \
-    .uvcalc_flag = UVCALC_TRANSFORM_CORRECT_SLIDE, \
-    .unwrapper = 1, \
+ \
+    .uvcalc_margin = _DNA_DEFAULT_ToolSettings_UVCalc_Margin, \
+    .uvcalc_flag = _DNA_DEFAULT_ToolSettings_UVCalc_Flag, \
+    .unwrapper = _DNA_DEFAULT_ToolSettings_UVCalc_Unwrapper, \
+    .uvcalc_iterations = _DNA_DEFAULT_ToolSettings_UVCalc_Iterations, \
+    .uvcalc_vertex_group = { 0 }, \
+    .uvcalc_vertex_group_factor = _DNA_DEFAULT_ToolSettings_UVCalc_VertexGroupFactor, \
+    .uvcalc_relative_scale = _DNA_DEFAULT_ToolSettings_UVCalc_RelativeScale, \
+    .uvcalc_allow_flips = _DNA_DEFAULT_ToolSettings_UVCalc_AllowFlips, \
+ \
     .select_thresh = 0.01f, \
  \
     .selectmode = SCE_SELECT_VERTEX, \

source/blender/makesdna/DNA_scene_types.h

@@ -1557,6 +1557,14 @@ typedef struct ToolSettings {
   char uv_sticky;
 
   float uvcalc_margin;
+  char uvcalc_allow_flips;
+  char _pad0[3];
+
+  int uvcalc_iterations;
+  float uvcalc_vertex_group_factor;
+  float uvcalc_relative_scale;
+
+  char uvcalc_vertex_group[64]; /* MAX_VGROUP_NAME */
 
   /* Auto-IK. */
   /** Runtime only. */
@@ -1584,7 +1592,7 @@ typedef struct ToolSettings {
   char gpencil_selectmode_edit;
   /** Stroke selection mode for Sculpt. */
   char gpencil_selectmode_sculpt;
-  char _pad0[6];
+  char _pad1[6];
 
   /** Grease Pencil Sculpt. */
   struct GP_Sculpt_Settings gp_sculpt;
@@ -1637,7 +1645,7 @@ typedef struct ToolSettings {
   short snap_flag_seq;
   short snap_flag_anim;
   short snap_uv_flag;
-  char _pad[4];
+  char _pad2[4];
   /** Default snap source, #eSnapSourceOP. */
   /**
    * TODO(@gfxcoder): Rename `snap_target` to `snap_source` to avoid previous ambiguity of
@@ -1701,7 +1709,7 @@ typedef struct ToolSettings {
 
   /** Normal Editing. */
   float normal_vector[3];
-  char _pad6[4];
+  char _pad3[4];
 
   /**
    * Custom Curve Profile for bevel tool:

source/blender/nodes/geometry/CMakeLists.txt

@@ -20,6 +20,9 @@ set(INC
   ../../modifiers
   ../../render
   ../../windowmanager
+  ../../../../extern/fmtlib/include
+  ../../../../intern/guardedalloc
+  ../../../../intern/slim
   # RNA_prototypes.h
   ${CMAKE_BINARY_DIR}/source/blender/makesrna
 )

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

@@ -67,12 +67,14 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
       mp_pin[i] = false;
       mp_select[i] = false;
     }
+
     geometry::uv_parametrizer_face_add(handle,
                                        face_index,
                                        face.size(),
                                        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

@@ -98,12 +98,14 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
       mp_pin[i] = false;
       mp_select[i] = false;
     }
+
     geometry::uv_parametrizer_face_add(handle,
                                        face_index,
                                        face.size(),
                                        mp_vkeys.data(),
                                        mp_co.data(),
                                        mp_uv.data(),
+                                       nullptr,
                                        mp_pin.data(),
                                        mp_select.data());
   });