均值坐标参数化（MVC Parameterization）

欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

均值坐标定义

均值坐标定义

$v_0是多边形v_1{v}_2{v}_3...v_n内的一点$

$就会存在均值坐标{\varphi }_i(v_0)=\frac{{\omega }_i}{{\sum }_{j=1}^n{\omega }_j}$

$其中{\omega }_i=\frac{tan\left(\frac{{\alpha }_{i-1}}{2}\right)+tan\left(\frac{{\alpha }_i}{2}\right)}{||v_i-v_0||}$

使得

$v_0=\sum _{i=1}^n{\varphi }_i(v_0)\cdot v_i\dots \dots （1）$

均值坐标参数化（MVC Parameterization）

该算法与tutte’s embedding 算法的基本过程一样，只是把线性组合的系数换成均值坐标表示。
该算法可以处理与圆盘同胚非闭合的三角网格面，他能保证3D与2D是一一映射的。

算法过程

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

行列式构建
$对（1）式两边都乘以\sum _{j\in \Omega (i)}{\omega }_{ij},再移项，得到$

$\left(\sum _{j\in \Omega (i)}{\omega }_{ij}\right)v_i-\sum _{j\in \Omega (i)}{\omega }_{ij}\cdot v_j=0$

利用上式，在具体实现的时候可以先计算出所有wij(没有连接关系则为0)，
然后对于边界点，行列式的值r[i][i]=1, b[i]=(ui, vi)
$对于非边界点r[i][i]=\left(\sum _{j\in \Omega (i)}{\omega }_{ij}\right)$

$r[i][j]=-{\omega }_{ij}，b[i]=(0,0)$

算法实现

代码链接点击前往
代码链接点击前往
代码链接点击前往

#include"include/PolyMesh/IOManager.h"
#include"include/PolyMesh/PolyMesh.h"
#include <Eigen/Sparse>
#define pi 3.1415926
using namespace acamcad;
using namespace polymesh;
using namespace std;
using namespace Eigen;

PolyMesh mesh;

void MVCCoordinates(PolyMesh* mesh)
{
	int v_n = mesh->numVertices();
	int iter;
	int boundary_num=0;

	VectorXd u(v_n), v(v_n);// u, v代表方程行列式的B列

	for (VertexIter v_it = mesh->vertices_begin(); v_it != mesh->vertices_end(); ++v_it)
	{
		if (mesh->isBoundary(*v_it))boundary_num++;
	}

	MHalfedge* first_heh= *mesh->halfedge_begin();
	// 查找边界
	for (HalfEdgeIter he_it = mesh->halfedge_begin(); he_it != mesh->halfedge_end(); ++he_it)
	{
		if (mesh->isBoundary(*he_it))
		{
			first_heh = *he_it;
			break;
		}
	}
	MHalfedge* iter_heh = first_heh->next();
	iter = 0;
	// 初始化边界点的uv值
	while (iter_heh != first_heh)
	{
		MVert* from_v = iter_heh->fromVertex();
		u[from_v->index()] = cos(double(2 * pi * iter / boundary_num));
		v[from_v->index()] = sin(double(2 * pi * iter / boundary_num));
		iter_heh = iter_heh->next();
		iter++;
	}
	u[first_heh->fromVertex()->index()] = cos(double(2 * pi * iter / boundary_num));
	v[first_heh->fromVertex()->index()] = sin(double(2 * pi * iter / boundary_num));

	SparseMatrix<double> weight(v_n, v_n);
	std::vector<Triplet<double>> triplet;

	// 针对每个面计算出边的wij. 后面可以聚合到矩阵中。
	for (FaceIter f_it = mesh->polyfaces_begin(); f_it != mesh->polyfaces_end(); ++f_it)
	{
		MHalfedge* heh = (*f_it)->halfEdge();
		MVert* v0 = heh->fromVertex();
		MVert* v1 = heh->toVertex();
		MHalfedge* next_heh = heh->next();
		MVert* v2 = next_heh->toVertex();

		double l2 = (v0->position() - v1->position()).norm();
		double l1 = (v0->position() - v2->position()).norm();
		double l0 = (v1->position() - v2->position()).norm();

		double angle = acos(dot(v2->position() - v0->position(), v1->position() - v0->position()) / (l1 * l2));
		triplet.push_back(Triplet<double>(v0->index(), v1->index(), tan(angle * 0.5) / l2));
		triplet.push_back(Triplet<double>(v0->index(), v2->index(), tan(angle * 0.5) / l1));

		angle = acos(dot(v0->position() - v1->position(), v2->position() - v1->position()) / (l2 * l0));
		triplet.push_back(Triplet<double>(v1->index(), v0->index(), tan(angle * 0.5) / l2));
		triplet.push_back(Triplet<double>(v1->index(), v2->index(), tan(angle * 0.5) / l0));

		angle = acos(dot(v0->position() - v2->position(), v1->position() - v2->position()) / (l0 * l1));
		triplet.push_back(Triplet<double>(v2->index(), v0->index(), tan(angle * 0.5) / l1));
		triplet.push_back(Triplet<double>(v2->index(), v1->index(), tan(angle * 0.5) / l0));
	}
	// 矩阵聚合
	weight.setFromTriplets(triplet.begin(), triplet.end());

	SparseMatrix<double> matrix(v_n, v_n);

	triplet.clear();
	for (VertexIter v_it = mesh->vertices_begin(); v_it != mesh->vertices_end(); ++v_it)
	{
		// 边界处只要设置 matirx[i][i]=1
		if (mesh->isBoundary(*v_it)) {
			triplet.push_back(Triplet<double>((*v_it)->index(), (*v_it)->index(), 1));
			continue;
		}

		// 内部点，设置 B列为0
		u[(*v_it)->index()] = 0;
		v[(*v_it)->index()] = 0;

		// maxtrix[i][i] = sum(wij), maxtrix[i][j] = -wij 
		for (VertexVertexIter vv_it = mesh->vv_iter(*v_it); vv_it.isValid(); ++vv_it)
		{
			triplet.push_back(Triplet<double>((*v_it)->index(), (*vv_it)->index(), -weight.coeff((*v_it)->index(), (*vv_it)->index())));
			triplet.push_back(Triplet<double>((*v_it)->index(), (*v_it)->index(), weight.coeff((*v_it)->index(), (*vv_it)->index())));
		}
	}
	matrix.setFromTriplets(triplet.begin(), triplet.end());
	SparseLU<SparseMatrix<double>> solver;
	solver.analyzePattern(matrix);
	solver.factorize(matrix);
	VectorXd result_u = solver.solve(u);
	VectorXd result_v = solver.solve(v);

	for (VertexIter v_it = mesh->vertices_begin(); v_it != mesh->vertices_end(); ++v_it)
	{
		(*v_it)->setPosition(result_u[(*v_it)->index()], result_v[(*v_it)->index()], 0);
	}
}

算法效果

原模型

参数化结果

原模型

参数化结果

---

本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。