网格参数化方法

参数化的作用

简单说就是把三维三角网格在二维上展开。由于曲率的存在，肯定存在变形。参数化要研究的问题就是如何让变形最小。

参数化后的二维平面可以存储很多信息，比如法向，纹理等。
可应用在纹理贴图，二维三维联动等。

tutte’s embedding 算法

该算法可以处理与圆盘同胚非闭合的三角网格面，他能保证3D与2D是一一映射的。

算法过程

1. 先把边界点按照顺序均匀地放置在圆边界上。
2. 对于内部的点，vi 是周围顶点的线性组合。
   $v_i=\sum _{j\in \Omega (i)}{d}_i{v}_j,d_i是顶点v_i的度。$
   n个未知顶点，n个方程刚好解出来。

算法实现

github: https://github.com/LightningBilly/ACMAlgorithms/tree/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw4

#include "hw4.h"
#include <fstream>
#include"../../PolyMesh/IOManager.h"
#include"../../PolyMesh/PolyMesh.h"
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <vector>
#include<iostream>
#include<math.h>
#include<algorithm>
// #include<omp.h>
#include <unistd.h>
using namespace acamcad;
using namespace polymesh;
using namespace std;


PolyMesh mesh;

#define M_PI 3.1415926

void tutte_parameterization() {
    char buffer[500];
    getcwd(buffer, 500);
    printf("The current directory is: %s/../\n", buffer);
    string mesh_path = buffer;
    mesh_path += "/../src/hw4/cow.obj";
    //读入网格
    // std::string mesh_path = argv[1];
    //PolyMesh *mesh1 = new PolyMesh();


    loadMesh(mesh_path, &mesh);

    //loadMesh("Bunny_head.obj", &mesh);
    mesh.updateMeshNormal();

    int F_N = mesh.numPolygons();
    int V_N = mesh.numVertices();

    // 计算出网格表面积，作为圆的面积
    double area_sum = 0;
    for (int i = 0; i < F_N; ++i)
    {
        auto f_h = mesh.polyface(i);
        auto itfv = mesh.fv_iter(f_h);
        auto v0 = (*itfv)->position();
        ++itfv;
        auto v1 = (*itfv)->position();
        ++itfv;
        auto v2 = (*itfv)->position();

        auto e0 = v1 - v0;
        auto e1 = v2 - v0;
        auto avec = cross(e0, e1);
        area_sum += avec.norm()/2.0;
    }

    // 将边界放置到圆上
    int boundary_num = 0;
    auto it1 = mesh.halfedge_begin();
    while (!mesh.isBoundary(*it1))
        it1++;
    auto he_start = *it1;
    auto he_it = he_start;
    do
    {
        he_it = (he_it)->next();
        boundary_num++;
    } while (he_it != he_start);

    double delta_angle = 2 * M_PI / boundary_num;
    double area_1_factor = sqrt(area_sum / M_PI);
    Eigen::VectorXd position_of_mesh; // 参数化坐标，前n个数字表示u[i], 后面n个数字表示v[i]
    position_of_mesh.resize(2 * V_N);
    for (int i = 0; i < boundary_num; ++i)
    {
        auto v_h = he_start->toVertex();
        position_of_mesh(v_h->index()) = area_1_factor * cos(i * delta_angle); // 设置u[i]
        position_of_mesh(v_h->index() + V_N) = area_1_factor * sin(-i * delta_angle); // 设置v[i]
        he_start = he_start->next();
    }


    // 填充Ax=b A矩阵和b微量
    typedef Eigen::Triplet<double> T;
    typedef Eigen::SparseMatrix<double> SMatrix;

    std::vector<T> tripletlist;
    Eigen::VectorXd bu = Eigen::VectorXd::Zero(V_N); // u向量
    Eigen::VectorXd bv = Eigen::VectorXd::Zero(V_N); // v向量
    // 遍历每个点，根据1领域情况填充
    for (auto it1 = mesh.vertices_begin(); it1 != mesh.vertices_end(); it1++)
    {
        int it1idx = (*it1)->index();
        if (mesh.isBoundary(*it1))
        {
            // 如果是边界说明已经有值了
            // 直接填充u[i]
            tripletlist.push_back(T(it1idx, it1idx, 1)); // 表示A[it1idx][it1idx]=1
            auto point = (*it1)->position();
            // 填充B列
            bu(it1idx) = position_of_mesh[it1idx];
            bv(it1idx) = position_of_mesh[it1idx + V_N];
        }
        else
        {
            // vi = sum(vj)/di, 移项后得 vi - sum(vj)/di=0, 都乘上度得 vi*di - sum(vj) = 0
            // 未知数b列为0
            for (auto it2 = mesh.vv_iter(*it1); it2.isValid() ; ++it2)
            {
                tripletlist.push_back(T(it1idx, (*it2)->index(), -1)); // A[i][j]=-1
            }
            tripletlist.push_back(T(it1idx, it1idx, mesh.valence(*it1))); // A[i][j] = di
        }
    }

    SMatrix coff(V_N, V_N);
    coff.setFromTriplets(tripletlist.begin(), tripletlist.end()); // 载入到矩阵
    Eigen::SparseLU<SMatrix> solver;
    solver.compute(coff);
    // 两列是独立的，可以分开求解
    Eigen::VectorXd xu = solver.solve(bu);
    Eigen::VectorXd xv = solver.solve(bv);

    Eigen::MatrixX2d uv;
    uv.resize(V_N, 2);

    //solve the positions
    uv.col(0) = xu;
    uv.col(1) = xv;

    for (int i = 0; i < V_N; i++)
    {
        auto v_h = mesh.vert(i);
        v_h->setPosition(uv(i, 0), uv(i, 1), 0);
    }
    writeMesh("bunny_tutte_para.obj", &mesh);
}

算法效果

最小二乘保角参数化（LSCM）

重心坐标

已知三角形的三点坐标，以及定义在顶点上的函数f，那么三角形内一点f值可以通过插值的方式得到。

$f(x)=\alpha f_i+\beta f_j+\gamma f_k$

其中

$\alpha =\frac{S_i}{S_i+S_j+S_k},S_i是x_i的对边与x组成的三角形的面积，S_j,S_k,\beta ,\gamma 同理$

线性分片性质

从上式可看出f(x)是一个分片线性函数。
那么, f(x)就可以由一个Jacobi矩阵乘上x再加一个列向量得到。
$f(x)=Jx+b，J是Jacobi矩阵$
$J就是先对f(x)求偏导，再对x求偏导$
$形式如下J=\left[\begin{array}{ccc}\frac{\partial u}{\partial x}& \frac{\partial v}{\partial x}& ...\\ \frac{\partial u}{\partial y}& \frac{\partial v}{\partial y}& ...\\ \frac{\partial u}{\partial z}& \frac{\partial v}{\partial z}& ...\end{array}\right]$
u, v 是f函数值的分量，如果有多个分量J就有多列。

梯度

上面的J矩阵，在三角形中表现为梯度。
梯度就是将f函数对x进行求导。
${\nabla }_{x}f(x)=f_i{\nabla }_x\alpha +f_j{\nabla }_x\beta +f_k\nabla \gamma$

面积Si要怎么求。
利用底乘以高。
底就是向量xk-xj的模，黄色向量是与xk-xj垂直的单位向量，高就是x-xj在黄向量的投影，可以用点积得到。
黄色向量可以使用xk-xj转90度取单位向量得到。
$\alpha =\frac{A_i}{A^T}=\frac{(x-x_j)\cdot \frac{(x_k-x_j{)}^{\perp }}{||x_k-x_j|{|}_2}||x_k-x_j|{|}_2}{2A^T}=\frac{(x-x_j)\cdot (x_k-x_j{)}^{\perp }}{2A^T}$

$只有前面和x变量有关，{\nabla }_x\alpha =\frac{(x_k-x_j{)}^{\perp }}{2A^T}$

${\nabla }_{x}f(x)=f_i\frac{(x_k-x_j{)}^{\perp }}{2A^T}+f_j\frac{(x_i-x_k{)}^{\perp }}{2A^T}+f_k\frac{(x_j-x_i{)}^{\perp }}{2A^T}$

这里可以看出一个三角的面的梯度是一个常数。

保角参数化基本思想

参数化以后是一个二维平面，保角的意思是希望二维平面的三角形与三维对应的三角形的角尽量相等。
那我们不妨先将三维的三角形摊平到2D上，然后进行变换。

相当于在2D三角形的顶点上定义了一个函数f, 经过函数f处理后就是参数化坐标，f(x) = [u, v]。

这里的x是一个二维坐标。
根据线性分片性质f(x) = Jx+b。

$J=\left[\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right]$

如果希望保角，那就要使J矩阵尽量是一个旋转矩阵。

$旋转矩阵是这样，\left[\begin{array}{cc}cos\partial & -sin\partial \\ sin\partial & cos\partial \end{array}\right]$
可以看出主对角是相等，副对角为相反数。

$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial y}，\frac{\partial u}{\partial y}=-\frac{\partial v}{\partial x}需要使这2个式子成立，完全相等是不可能的，只能是尽量接近$

可以建立能量函数$Q=\sum Q_t=\sum S_t\left[{\left(\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2\right]$
一个t代表一个三角形。

优化这个函数使得Q最小就可以保角了。

下面求解一下J矩阵。

$J={\nabla }_{x}f(x)=f_i\frac{(x_k-x_j{)}^{\perp }}{2A^T}+f_j\frac{(x_i-x_k{)}^{\perp }}{2A^T}+f_k\frac{(x_j-x_i{)}^{\perp }}{2A^T}$
先解决旋转问题。
$旋转90度矩阵：rotate90=\left[\begin{array}{cc}cos90^0& -sin90^0\\ sin90^0& cos90^0\end{array}\right]=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$

$(x_k-x_j{)}^{\perp }=\left[\begin{array}{cc}cos90^0& -sin90^0\\ sin90^0& cos90^0\end{array}\right]\left[\begin{array}{c}{x}_k-x_j\\ y_k-y_j\end{array}\right]=\left[\begin{array}{c}-(y_k-y_j)\\ x_k-x_j\end{array}\right]$

其他同理。

J矩阵坐标表示：

$J=\frac{1}{2A}\left[\begin{array}{ccc}-(y_k-y_j)& -(y_i-y_k)& -(y_j-y_i)\\ x_k-x_j& x_i-x_k& x_j-y_i\end{array}\right]\left[\begin{array}{c}{f}_i\\ f_j\\ f_k\end{array}\right]$

$=\frac{1}{2A}\left[\begin{array}{ccc}-(y_k-y_j)& -(y_i-y_k)& -(y_j-y_i)\\ x_k-x_j& x_i-x_k& x_j-y_i\end{array}\right]\left[\begin{array}{cc}{u}_i& v_i\\ u_j& v_j\\ u_k& v_k\end{array}\right]=\left[\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right]$

$=\frac{1}{2A}\left[\begin{array}{ccc}{k}_1& k_2& k_3\\ k_4& k_5& k_6\end{array}\right]\left[\begin{array}{cc}{u}_i& v_i\\ u_j& v_j\\ u_k& v_k\end{array}\right]$

$=\left[\begin{array}{ccc}{k}_1& k_2& k_3\\ k_4& k_5& k_6\end{array}\right]\left[\begin{array}{cc}{u}_i& v_i\\ u_j& v_j\\ u_k& v_k\end{array}\right]$

$=\left[\begin{array}{cc}{k}_1\ast u_i+k_2\ast u_j+k_3\ast u_k& k_1\ast v_i+k_2\ast v_j+k_3\ast v_k\\ k_4\ast u_i+k_5\ast u_j+k_6\ast u_k& k_4\ast v_i+k_5\ast v_j+k_6\ast v_k\end{array}\right]$

$$\begin{gathered} k_n是已知量，代入Q式子中，Q是关于u,v的二次多项， \\ 只要求驻点就可以得到最小值。求导后可以通过最小二乘解决。 \end{gathered}$$

代码实现

实际实现时是用了复分析的方法，具体可以参考以下文章：
https://blog.csdn.net/why18767183086/article/details/107870019
https://blog.csdn.net/weixin_49732532/article/details/107822729
https://www.zhihu.com/people/allan-35-49/posts
https://github.com/icemiliang/lscm

效果

ABF

Angle-Based Flattening(ABF) 是计算出参数化后各三角形的角度，使得角度与原来的角度相差最小。然后再利用角度来恢复出三角形的参数化方法。

基于这个思想给出能量方程以及约束条件：

$$E_{ABF}=\sum _t\sum _{i=1}^3{\omega }_i^t({\alpha }_i^t-{\beta }_i^t{)}^2$$

$$\begin{gathered} {\beta }_i^t是优化的目标角度，{\alpha }_i^t是变量 \\ {\omega }_i^t=({\beta }_i^t{)}^{-2} \end{gathered}$$

$${\beta }_i^t需要归一化,{\beta }_i^t=\{\begin{array}{l}\frac{{\tilde{\beta }}_i^t\cdot 2\pi }{{\sum }_i{\beta }_i^t},内部结点\\ {\tilde{\beta }}_i^t,边界结点\end{array}$$

约束条件如下：
所有的角都要为正
${\alpha }_i^t>0$

三角形内角和为180度
${\alpha }_1^t+{\alpha }_2^t+{\alpha }_3^t=\pi$

一个项点所有的所有三角形，以该点所形成的角和为360度
$\sum _{t\in \Omega (v)}{\alpha }_k^t=2\pi$

限制边长关系使得最后三角形可以合上
${\prod }_{t\in \Omega (v)}\sin {\alpha }_{k\oplus 1}^t={\prod }_{t\in \Omega (v)}\sin {\alpha }_{k⊝1}^t$

LABF

上述方法中
${\prod }_{t\in \Omega (v)}\sin {\alpha }_{k\oplus 1}^t={\prod }_{t\in \Omega (v)}\sin {\alpha }_{k⊝1}^t$
这个条件是一个非线性约束，这使得求解时非常麻烦。
LABF的思想是先把这个条件变成一个线性约束条件，再来求解。

$设{\alpha }_i^t={\beta }_i^t+e_i^t$

先对原式取log
$$log\left(\prod _{t\in \Omega (v)}\sin {\alpha }_{k\oplus 1}^t\right)=log\left(\prod _{t\in \Omega (v)}\sin {\alpha }_{k⊝1}^t\right)$$

$$\sum _{t\in \Omega (v)}\log (\sin {\alpha }_{k\oplus 1}^t)=\sum _{t\in \Omega (v)}\log (\sin {\alpha }_{k⊝1}^t)$$

利用泰勒展开得到

$$\log (\sin {\alpha }_{k\oplus 1}^t)=\log (\sin {\beta }_{k\oplus 1}^t+e_{k\oplus 1}^t)$$

$$=\log (\sin {\beta }_{k\oplus 1}^t)+e_{k\oplus 1}^t\cot {\beta }_{k\oplus 1}^t+...$$

那么原来的E可以转化成 e的方程
$$E_{ABF}=\sum _t\sum _{i=1}^3{\omega }_i^t(e_i^t{)}^2$$

利用拉格朗日乘子法可以得到

$$Ae=b$$
$$\left[\begin{array}{cc}D& A^T\\ A& 0\end{array}\right]\left[\begin{array}{c}e\\ \lambda \end{array}\right]=\left[\begin{array}{c}0\\ b\end{array}\right]$$

解得：

$$e=D^{-1}{A}^T(A{D}^{-1}{A}^T{)}^{-1}b$$

基于角度还原出参数化结果

贪心算法
先放一个种子三角形，对3条边进行深度优先遍历摆放三角形。
该方法不足之处是，由于浮点有计算误差，摆放过程中误差会越来越大，最后导致有一些三角形的边无法合上。

基于优化的算法

根据2边之比等于对角的sin之比可以得到
$\frac{||\overrightarrow{P_1{P}_3}||}{||\overrightarrow{P_1{P}_2}||}=\frac{\sin {\alpha }_2}{\sin {\alpha }_3}$

$向量\overrightarrow{P_1{P}_3}可以看成是\overrightarrow{P_1{P}_2}逆时针旋转再乘以一个比例得到$

$\overrightarrow{P_1{P}_3}=\frac{\sin {\alpha }_2}{\sin {\alpha }_3}\left[\begin{array}{cc}\cos {\alpha }_1& -\sin {\alpha }_1\\ \sin {\alpha }_1& \cos {\alpha }_1\end{array}\right]\overrightarrow{P_1{P}_2}$

$所以对于一个正常的三角形\overrightarrow{P_1{P}_3}=\frac{\sin {\alpha }_2}{\sin {\alpha }_3}\left[\begin{array}{cc}\cos {\alpha }_1& -\sin {\alpha }_1\\ \sin {\alpha }_1& \cos {\alpha }_1\end{array}\right](P_2-P_1)-(P_3-P_1)=0$

实际中做不到等于0，只能是使其越接近越好。

可以针对所有三角形列出一个能量方程

$$E=\sum _t||\frac{\sin {\alpha }_2^t}{\sin {\alpha }_3^t}\left[\begin{array}{cc}\cos {\alpha }_1^t& -\sin {\alpha }_1^t\\ \sin {\alpha }_1^t& \cos {\alpha }_1^t\end{array}\right](P_2^t-P_1^t)-(P_3^t-P_1^t)|{|}^2$$

优化上述方程即可得到参数化结果。

as-rigid-as-possible (ASAP)

尽量做保边长的变换。就是为每个三个角形找到一个单纯的旋转矩阵，lscm是找到一个旋转+缩放。asap只要旋转，不要缩放。

能量方程：
$$E(u,L)=\sum _t{A}_t||J_t-L_t|{|}_F^2$$

$L_t是三角形的目标旋转矩阵$
$J_t是三角形的实际的变换矩阵$

这里有2个变量可以使用交替迭代的方法求解。

$当给定J_t时，由于每个三角形是独立的，所以L_t可分开求解$
$对J_t做SVD分解，J_t=U\sum V^T,\sum =\left[\begin{array}{cc}{\sigma }_1& 0\\ 0& {\sigma }_2\end{array}\right],L_t最优值为U{V}^T$

$当L_t确定时，E(u,L)是一个二次多项式，可通过求导，利用最小二乘求解。$

具体求解时E可以表示成边和点的多项式。可以参考https://blog.csdn.net/why18767183086/article/details/108034725

代码实现

// arap_parameterizations
#include <fstream>
#include"../../PolyMesh/IOManager.h"
#include"../../PolyMesh/PolyMesh.h"
#include<Eigen/SVD>
#include<Eigen/Dense>
#include<Eigen/Sparse>
#include <vector>
#include<iostream>
#include<math.h>
#include<algorithm>
#include <unistd.h>
#include "hw5.h"
using namespace acamcad;
using namespace polymesh;
using namespace std;

PolyMesh mesh;

void Tutte(Eigen::MatrixX2d& uv)
{
    int F_N = mesh.numPolygons();
    int V_N = mesh.numVertices();

    //calc surface area
    double area_sum = 0;
    for (int i = 0; i < F_N; ++i)
    {
        auto f_h = mesh.polyface(i);
        auto itfv = mesh.fv_iter(f_h);
        auto v0 = (*itfv)->position();
        ++itfv;
        auto v1 = (*itfv)->position();
        ++itfv;
        auto v2 = (*itfv)->position();

        auto e0 = v1 - v0;
        auto e1 = v2 - v0;
        auto avec = cross(e0, e1);
        area_sum += avec.norm()/2.0;
    }

    //set the boundary vertices to circle
    int boundary_num = 0;
    auto it1 = mesh.halfedge_begin();
    while (!mesh.isBoundary(*it1))
        it1++;
    auto he_start = *it1;
    auto he_it = he_start;
    do
    {
        he_it = (he_it)->next();
        boundary_num++;
    } while (he_it != he_start);

    double delta_angle = 2 * M_PI / boundary_num;
    double area_1_factor = sqrt(area_sum / M_PI);
    Eigen::VectorXd position_of_mesh;
    position_of_mesh.resize(2 * V_N);
    for (int i = 0; i < boundary_num; ++i)
    {
        auto v_h = he_start->toVertex();
        position_of_mesh(v_h->index()) = area_1_factor * cos(i * delta_angle);
        position_of_mesh(v_h->index() + V_N) = area_1_factor * sin(-i * delta_angle);
        he_start = he_start->next();
    }


    //calc the matrix
    typedef Eigen::Triplet<double> T;
    typedef Eigen::SparseMatrix<double> SMatrix;

    std::vector<T> tripletlist;
    Eigen::VectorXd bu = Eigen::VectorXd::Zero(V_N);
    Eigen::VectorXd bv = Eigen::VectorXd::Zero(V_N);
    for (auto it1 = mesh.vertices_begin(); it1 != mesh.vertices_end(); it1++)
    {
        int it1idx = (*it1)->index();
        if (mesh.isBoundary(*it1))
        {
            tripletlist.push_back(T(it1idx, it1idx, 1));
            auto point = (*it1)->position();
            bu(it1idx) = position_of_mesh[it1idx];
            bv(it1idx) = position_of_mesh[it1idx + V_N];
        }
        else
        {
            for (auto it2 = mesh.vv_iter(*it1); it2.isValid() ; ++it2)
            {
                tripletlist.push_back(T(it1idx, (*it2)->index(), -1));
            }
            tripletlist.push_back(T(it1idx, it1idx, mesh.valence(*it1)));
        }
    }

    SMatrix coff(V_N, V_N);
    coff.setFromTriplets(tripletlist.begin(), tripletlist.end());
    Eigen::SparseLU<SMatrix> solver;
    solver.compute(coff);
    Eigen::VectorXd xu = solver.solve(bu);
    Eigen::VectorXd xv = solver.solve(bv);

    //solve the inner positions
    uv.col(0) = xu;
    uv.col(1) = xv;
}


void arap_parameterization()
{
    /*
    if (argc != 3)
    {
        std::cout << "========== Hw5 Usage  ==========\n";
        std::cout << std::endl;
        std::cout << "Input:	ACAM_mesh_HW5.exe	input_mesh.obj	output_mesh.obj\n";
        std::cout << std::endl;
        std::cout << "=================================================\n";
        return;
    }
    */
    char buffer[500];
    getcwd(buffer, 500);
    printf("The current directory is: %s/../\n", buffer);
    string mesh_path = buffer;
    mesh_path += "/../src/hw5/cow.obj";
    //读入网格
    loadMesh(mesh_path, &mesh);

    std::string out_path = buffer;
    out_path += "/../src/hw5/output_mesh.obj";

    //loadMesh("cow.obj", &mesh);
    mesh.updateMeshNormal();
    int nf = mesh.numPolygons();
    int nv = mesh.numVertices();

    //calc local coordinate
    Eigen::MatrixXd localcoord;
    localcoord.resize(nf, 6);

#pragma omp parallel for
    for (int i = 0; i < nf; i++)
    {
        auto f_h = mesh.polyface(i);

        auto n = f_h->normal();
        auto itfv = mesh.fv_iter(f_h);
        auto v0 = (*itfv)->position();
        ++itfv;
        auto v1 = (*itfv)->position();
        ++itfv;
        auto v2 = (*itfv)->position();

        MVector3 e = v1 - v0;

        auto x_ = e.normalized();
        auto y_ = cross(n,x_);

        auto e1 = v2 - v0;

        localcoord.row(i) << 0, 0, e.norm(), 0, dot(e1 , x_), dot(e1 , y_);
    }


    //calc cot-weight laplacian matrix
    vector<double> cots(mesh.numHalfEdges());
    vector<Eigen::Triplet<double>> trivec;
    for (auto ithe = mesh.halfedge_begin(); ithe != mesh.halfedge_end(); ithe++)
    {
        if (mesh.isBoundary(*ithe))
            continue;

        auto v1 = (*ithe)->fromVertex()->position();
        auto v2 = (*ithe)->toVertex()->position();
        auto v0 = (*ithe)->next()->toVertex()->position();

        auto e0 = v1 - v0;
        auto e1 = v2 - v0;
        double cotangle = dot(e0, e1) / cross(e0, e1).norm();

        cots[(*ithe)->index()] = cotangle;
        int vid0 = (*ithe)->fromVertex()->index();
        int vid1 = (*ithe)->toVertex()->index();
        trivec.emplace_back(vid0, vid0, cotangle);
        trivec.emplace_back(vid1, vid1, cotangle);
        trivec.emplace_back(vid0, vid1, -cotangle);
        trivec.emplace_back(vid1, vid0, -cotangle);
    }
    Eigen::SparseMatrix<double> smat;
    smat.resize(nv, nv);
    smat.setFromTriplets(trivec.begin(), trivec.end());

    Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
    solver.compute(smat);


    Eigen::MatrixX2d uv;
    uv.resize(nv, 2);

    //tutte initialization;
    Tutte(uv);

    vector<Eigen::Matrix2d> Lts;
    Lts.resize(nf);
    //local-global iteration
    for (int iter = 0; iter < 100; iter++)
    {
        //local calc Lt
#pragma omp parallel for
        for (int i = 0; i < nf; i++)
        {
            auto f_h = mesh.polyface(i);
            auto itfv = mesh.fv_iter(f_h);
            auto i0 = (*itfv)->index();
            ++itfv;
            auto i1 = (*itfv)->index();
            ++itfv;
            auto i2 = (*itfv)->index();

            Eigen::Matrix2d P, S, J;

            P << uv(i1, 0) - uv(i0, 0), uv(i2, 0) - uv(i0, 0), uv(i1, 1) - uv(i0, 1), uv(i2, 1) - uv(i0, 1);
            S << localcoord(i, 2) - localcoord(i, 0), localcoord(i, 4) - localcoord(i, 0),
                    localcoord(i, 3) - localcoord(i, 1), localcoord(i, 5) - localcoord(i, 1);

            J = P * S.inverse();

            Eigen::JacobiSVD<Eigen::Matrix2d> svd(J, Eigen::ComputeFullU | Eigen::ComputeFullV);

            Eigen::Matrix2d U = svd.matrixU();
            Eigen::Matrix2d V = svd.matrixV();

            Eigen::Matrix2d R = U * V.transpose();

            if (R.determinant() < 0)
            {
                U(0, 1) = -U(0, 1);
                U(1, 1) = -U(1, 1);
                R = U * V.transpose();
            }
            Lts[i] = R;
        }

        //global calc b
        Eigen::VectorXd bu, bv;
        bu.setZero(nv);
        bv.setZero(nv);

        for (int i = 0; i < nf; i++)
        {
            auto f_h = mesh.polyface(i);

            auto n = f_h->normal();
            auto itfv = mesh.fv_iter(f_h);

            auto i0 = (*itfv)->index();
            ++itfv;
            auto i1 = (*itfv)->index();
            ++itfv;
            auto i2 = (*itfv)->index();

            auto he2 = f_h->halfEdge();
            auto he0 = he2->next();
            auto he1 = he0->next();

            Eigen::Vector2d e0, e1, e2;
            e0 << localcoord(i, 2), localcoord(i, 3);
            e1 << localcoord(i, 4) - localcoord(i, 2), localcoord(i, 5) - localcoord(i, 3);
            e2 << -localcoord(i, 4), -localcoord(i, 5);
            Eigen::Vector2d b0 = cots[he0->index()] * Lts[i] * e0;
            bu[i0] -= b0[0];
            bv[i0] -= b0[1];

            bu[i1] += b0[0];
            bv[i1] += b0[1];
            Eigen::Vector2d b1 = cots[he1->index()] * Lts[i] * e1;
            bu[i1] -= b1[0];
            bv[i1] -= b1[1];

            bu[i2] += b1[0];
            bv[i2] += b1[1];

            Eigen::Vector2d b2 = cots[he2->index()] * Lts[i] * e2;
            bu[i2] -= b2[0];
            bv[i2] -= b2[1];

            bu[i0] += b2[0];
            bv[i0] += b2[1];

        }

        //global solve
        uv.col(0) = solver.solve(bu);
        uv.col(1) = solver.solve(bv);
    }

    for (int i = 0; i < nv; i++)
    {
        auto v_h = mesh.vert(i);
        v_h->setPosition(uv(i, 0), uv(i, 1), 0);
    }

    writeMesh(out_path, &mesh);

}

算法效果

原型图

参数化图