【数字几何处理】参数化：Tutte‘s embeddingLeast Squares Conformal Mappings 源码+介绍

Tutte’s embedding

原理：如果边界顶点有序的落在一个凸多边形上，且内部的顶点是其邻居的线性组合，那么（u，v）坐标参数化是双射的。

$[1000⋯0100⋯⋯⋯⋯⋯⋯wi,0wi,1−Σwi,jwi,4⋯⋯⋯⋯⋯⋯][u0v0u1v1u2v2u3v3u4v4]=[u0Boundv0Boundu1Boundv1Bound000000]\begin{bmatrix} 1&0&0&0&\cdots \\ 0&1&0&0&\cdots\\ \cdots & \cdots & \cdots & \cdots & \cdots\\ w_{i,0}&w_{i,1}&-\Sigma w_{i,j}&w_{i,4} & \cdots\\ \cdots & \cdots & \cdots & \cdots & \cdots \end{bmatrix}\begin{bmatrix}u0&v0\\ u1&v1\\u2&v2\\u3&v3\\u4&v4\\\end{bmatrix}=\begin{bmatrix}u0Bound&v0Bound\\ u1Bound&v1Bound\\0&0\\0&0\\0&0\\\end{bmatrix}$

对于边界上的点，uv的值已知，对于内部的点，其和邻居的线性组合为0。
当边ij存在时， $w_{ij}$ 不为0，反之， $w_{ij}$ 为0。

$w_{ij}$ 有许多种取值方法：
$1.uniformwij=11.uniform\space w_{ij}=1$
$2.edge−basedwij=1∣∣vi−vj∣∣2.edge-based\space w_{ij}=\frac{1}{||v_{i}-v_{j}||}$
$3.harmonicwij=cotαij+βij23.harmonic\space w_{ij}=\frac{cot\alpha_{ij}+\beta_{ij}}{2}$
$4.mean−valuewij=tan(γij)+tan(δij)2∣∣vi−vj∣∣4.mean-value\space w_{ij}=\frac{tan(\gamma_{ij})+tan(\delta_{ij})}{2||v_{i}-v_{j}||}$

代码步骤：
1.找出边界
2.把边界映射到二维的圆上
3.计算矩阵L
3.1若该点为边界点，跳过
3.2若该点不是边界点，计算权重
3.3将所有的边界点 $L (i, i)$ 的值设为0
4.求解方程组
4.1求解 $Lu=b_{u}$
4.2求解 $Lv=b_{v}$

对于步骤1：
1.要考虑有多个边界的情况：
1.1找出所有边界点
1.2将边界点归类到各自的边界中
1.3以最长的边界作为结果

#include "tutte.h"
#include "map_vertices_to_circle.h"
#include "boundary_loop.h"
#include "cotmatrix.h"void tutte(const Eigen::MatrixXd & V,const Eigen::MatrixXi & F,Eigen::MatrixXd & U)
{Eigen::VectorXi boundary;cgpl::boundary_loop(F, boundary);Eigen::MatrixXd UV;cgpl::map_vertices_to_circle(V, boundary, UV);std::vector<bool> isBound(V.rows(), false);for (int i = 0; i < boundary.rows(); i++){isBound[boundary(i, 0)] = true;}Eigen::SparseMatrix<double> L;typedef Eigen::Triplet<double> T;std::vector<T> tripletList;for (int f = 0; f < F.rows(); f++) {for (int i = 0; i < 3; i++) {int ii = F(f, i);if (isBound[ii]){continue;}int j = (i + 1) % 3;int k = (i + 2) % 3;double norm_ij = (V.row(F(f, i)) - V.row(F(f, j))).norm();double norm_ik = (V.row(F(f, i)) - V.row(F(f, k))).norm();// case i != jtripletList.push_back(T(F(f, i), F(f, j), 1.0/norm_ij));tripletList.push_back(T(F(f, i), F(f, k), 1.0/norm_ik));// case i == jtripletList.push_back(T(F(f, i), F(f, i), -1.0/norm_ij));tripletList.push_back(T(F(f, i), F(f, i), -1.0/norm_ik));}}for (int i = 0; i < boundary.rows(); i++){tripletList.push_back(T(boundary(i), boundary(i), 1));}L.resize(V.rows(), V.rows());L.setFromTriplets(tripletList.begin(), tripletList.end());U.resize(V.rows(), 2);Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>solver;Eigen::VectorXd Bu = Eigen::VectorXd::Zero(V.rows());Eigen::VectorXd Bv = Eigen::VectorXd::Zero(V.rows());for (int i = 0; i < boundary.rows(); i++){Bu(boundary[i]) = UV(i, 0);Bv(boundary[i]) = UV(i, 1);}solver.compute(L);Eigen::VectorXd u=solver.solve(Bu);Eigen::VectorXd v=solver.solve(Bv);U << u, v;
}

#include "boundary_loop.h"
#include "triangle_triangle_adjacency.h"
#include "vertex_triangle_adjacency.h"
#include "is_border_vertex.h"
#include <set>
void cgpl::boundary_loop(const Eigen::MatrixXi & F,std::vector<std::vector<int> >& L)
{using namespace std;using namespace Eigen;if(F.rows() == 0)return;MatrixXi TT;vector<std::vector<int>> VF, VFi;triangle_triangle_adjacency(F,TT);vertex_triangle_adjacency(F.maxCoeff()+1,F,VF,VFi);vector<bool> unvisited = is_border_vertex(F);set<int> unseen;for (size_t i = 0; i < unvisited.size(); ++i){if (unvisited[i])unseen.insert(unseen.end(),i);}while (!unseen.empty()){vector<int> l;// Get first vertex of loopint start = *unseen.begin();unseen.erase(unseen.begin());unvisited[start] = false;l.push_back(start);bool done = false;while (!done){// Find next vertexbool newBndEdge = false;int v = l[l.size()-1];int next;for (int i = 0; i < (int)VF[v].size() && !newBndEdge; i++){int fid = VF[v][i];if (TT.row(fid).minCoeff() < 0.) // Face contains boundary edge{int vLoc = -1;if (F(fid,0) == v) vLoc = 0;if (F(fid,1) == v) vLoc = 1;if (F(fid,2) == v) vLoc = 2;int vNext = F(fid,(vLoc + 1) % F.cols());newBndEdge = false;if (unvisited[vNext] && TT(fid,vLoc) < 0){next = vNext;newBndEdge = true;}}}if (newBndEdge){l.push_back(next);unseen.erase(next);unvisited[next] = false;}elsedone = true;}L.push_back(l);}}void cgpl::boundary_loop(const Eigen::MatrixXi & F,std::vector<int> & L)
{using namespace Eigen;using namespace std;if(F.rows() == 0)return;vector<vector<int> > Lall;boundary_loop(F,Lall);int idxMax = -1;size_t maxLen = 0;for (size_t i = 0; i < Lall.size(); ++i){if (Lall[i].size() > maxLen){maxLen = Lall[i].size();idxMax = i;}}//Check for meshes without boundaryif (idxMax == -1){L.clear();return;}L.resize(Lall[idxMax].size());for (size_t i = 0; i < Lall[idxMax].size(); ++i){L[i] = Lall[idxMax][i];}
}void cgpl::boundary_loop(const Eigen::MatrixXi & F,Eigen::VectorXi & L)
{using namespace Eigen;using namespace std;if(F.rows() == 0)return;vector<int> Lvec;boundary_loop(F,Lvec);L.resize(Lvec.size());for (size_t i = 0; i < Lvec.size(); ++i)L(i) = Lvec[i];
}

Least Squares Conformal Mappings

原理：尽量保角和保面积

保角

$Jt=[∂u∂x∂u∂y∂v∂x∂v∂y]=[a−bba]=s[cosθ−sinθsinθcosθ]J_{t}=\begin{bmatrix} \frac{\partial u}{\partial x}&\frac{\partial u}{\partial y}\\\frac{\partial v}{\partial x}&\frac{\partial v}{\partial y}\end{bmatrix}=\begin{bmatrix}a&-b\\b&a\end{bmatrix}=s \begin{bmatrix}cos\theta&-sin\theta\\sin\theta&cos\theta\end{bmatrix}$

得出：

$∂u∂x=∂v∂y∂u∂y=−∂v∂x\frac{\partial u}{\partial x}=\frac{\partial v}{\partial y}\space \frac{\partial u}{\partial y}=-\frac{\partial v}{\partial x}$

$▽u=▽v⊥\triangledown u=\triangledown v^{\perp}$

即：

$minu,v12∫S∣∣▽u−▽v⊥∣∣2dAmin_{u,v}\frac{1}{2}\int_{S}||\triangledown u-\triangledown v^{\perp}||^{2}dA$

$minu,v∫S(12∣∣▽u∣∣2+12∣∣▽v∣∣2−▽u⋅▽v⊥)dAmin_{u,v}\int_{S}(\frac{1}{2}||\triangledown u||^{2}+\frac{1}{2}||\triangledown v||^{2}-\triangledown u \cdot \triangledown v^{\perp})dA$

$minUUT([LL]−A)Umin_{U}U^{T}(\begin{bmatrix}L&\\ &L\end{bmatrix}-A)U$

L是laplacian matrix，A是vector area，把L放在对角线上（专业术语叫啥忘记了）可以只用一次就可计算出u和v的值。

$detJ_{t}=1$

$UT([MM])U=1U^{T}(\begin{bmatrix}M&\\ &M\end{bmatrix})U=1$

求解

$minv12vTAvsubjecttovTBv=1min_{v}\frac{1}{2} v^{T} Av \space subject\space to\space v^{T} Bv= 1$

用拉格朗日乘数法把它转化成广义特征值问题

$12vTAv+λ(1−vTBv)\frac{1}{2} v^{T} Av +\lambda(1- v^{T} Bv)$

$Av−λBv=0Av-\lambda Bv=0$

$Av=λBvAv=\lambda Bv$

$1- v^{T} Bv=0$

$v^{T} Bv=1$

结果中前两个为平凡解，所以我们取第三个作为结果

#include "lscm.h"
#include "vector_area_matrix.h"#include <Eigen/SVD>
#include "eigs.h"
#include "repdiag.h"
#include "cotmatrix.h"
#include "massmatrix.h"
void lscm(const Eigen::MatrixXd & V,const Eigen::MatrixXi & F,Eigen::MatrixXd & U)
{// Solve optimization as a generalized Eigen value problem//      min_U U'QU  subject to  U'BU = 1int n = V.rows(); // Compute Q & BEigen::SparseMatrix<double> A, L, Q;cgpl::cotmatrix(V, F, L);cgpl::repdiag(L, 2, Q);vector_area_matrix(F, A);Q = Q - A;Eigen::SparseMatrix<double> M, B;cgpl::massmatrix(V, F, M);cgpl::repdiag(M, 2, B);// SolveEigen::MatrixXd sU;Eigen::VectorXd sS;cgpl::eigs(Q, B, 4, sU, sS);// Somehow, first 2 eigen value ~e^{-13}//      3, 4, ... looks more reasonableU.resize(n, 2);U << sU.col(2).topRows(n), sU.col(2).bottomRows(n);
}

结果

Tutte’s embedding

Least Squares Conformal Mappings