三角网格的拉普拉斯矩阵(Laplacian Matrix)的并行计算

这篇文章基于OpenMesh、Eigen和OpenMP，来介绍如何并行地计算三角网格常见的拉普拉斯矩阵(Laplacian Matrix)

在网格参数化中，拉普拉斯矩阵是一个出现频率很高的字眼，相当一部分的参数化方法都直接或间接地与拉普拉斯矩阵有关，就我实际使用的就有：

图特映射(Tutte Embedding)
调和映射(Harmonic Mapping)
里奇流(Ricci Flow)
卡拉比流(Calabi Flow)
李流(Lie Advection)
……

这篇文章基于OpenMesh、Eigen和OpenMP，来介绍如何并行地计算三角网格的拉普拉斯矩阵(Laplacian Matrix)。

拉普拉斯在三角网格上有很多种形式，我们将介绍最基本的形式和最常见的形式。

简单拉普拉斯

如果只考虑三角网格的拓扑关系，那么三角网格的拉普拉斯矩阵和简单图的拉普拉斯是一样的。参照维基百科，我们有

$L=D-A$

即[拉普拉斯矩阵]=[度矩阵]-[邻接矩阵]。另一个我喜欢的等效的定义是：

$L_{ij}=\left\{\begin{matrix} -1 &, e_{ij} \in E \\ \sum_{e_{ik} \in E}^{ }-L_{ik} & , i=j \\ 0 & , otherwise \end{matrix}\right.$

这个定义比较易于计算。可以看出，拉普拉斯矩阵一般是一个稀疏的正定矩阵，在Eigen中适合使用SparseMatrix存储， $L_{ij}$ 当且仅当 $i=j$ 或 $e_{ij} \in E$ 时非0。下面的C++例子展示了如何从一个OpenMesh网格创建对应的基本形式的拉普拉斯矩阵：

//在SparseMatrix中线程安全地为Laplacian分配空间
template<int IsRowMajor>
Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor> PrepareLaplacian(const MyMesh & mesh)
{
	const int vertices = mesh.n_vertices();
	Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor> A(vertices, vertices);

	//遍历网格点-点连接，同时：
	//1、缓存每个点的邻居点
	//2、为SparseMatrix计算需要预分配的空间
	VectorXi sizes(vertices);
	std::vector<std::vector<MyMesh::VertexHandle>> neighbors(vertices);
#pragma omp parallel for schedule(dynamic, 1000)
	for (int i = 0; i < vertices; i++) {
		neighbors[i].reserve(10);
		for (auto vv : mesh.vv_range(MyMesh::VertexHandle(i))) {
			neighbors[i].push_back(vv);
		}
		sizes[i] = neighbors[i].size() + 1;		//多出一个对角线元素
	}

	//为行列分配元素空间
	A.reserve(sizes);
#pragma omp parallel for schedule(dynamic, 1000)
	for (int i = 0; i < vertices; i++) {
		for (int j = 0; j < neighbors[i].size(); j++) {
			if (A.IsRowMajor) {
				A.insert(i, neighbors[i][j].idx());
			}
			else {
				A.insert(neighbors[i][j].idx(), i);
			}
		}
		A.insert(i, i);
	}

	return A;
}

//根据非对角线元素，更新对角线元素
template<int IsRowMajor>
void UpdateDiagonal(Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor>& A)
{
#pragma omp parallel for schedule(dynamic, 1000)
	for (int i = 0; i < A.outerSize(); i++)
	{
		A.coeffRef(i, i) = 0;
		MyMesh::Scalar sum = 0.0;
		for (Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor>::InnerIterator it(A, i); it; ++it)
		{
			sum += it.value();
		}
		A.coeffRef(i, i) = -sum;
	}
}

//更新基本形式的拉普拉斯矩阵
template<int IsRowMajor>
void UpdateTutte(const MyMesh & mesh, Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor>& A)
{
	//计算非对角线元素
#pragma omp parallel for schedule(dynamic, 100)
	for (int i = 0; i < mesh.n_edges(); i++) {
		//      v1
		//    //||\\
		//   v2 || v3    
		//    \\||//
		//      v0
		const auto eh = mesh.edge_handle(i);
		const auto he01 = mesh.halfedge_handle(eh, 0);
		const auto v0 = mesh.from_vertex_handle(he01);
		const auto v1 = mesh.to_vertex_handle(he01);

		A.coeffRef(v0.idx(), v1.idx()) = -1;
		A.coeffRef(v1.idx(), v0.idx()) = -1;
	}

	//计算对角线元素
	UpdateDiagonal(A);
}

//构建基本形式的拉普拉斯矩阵
template<int IsRowMajor>
Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor> BuildTutte(const MyMesh & mesh)
{
	auto A = PrepareLaplacian<IsRowMajor>(mesh);
	UpdateTutte(mesh, A);
	return A;
}

函数比较多，主要是为了实际应用时的代码复用。最后的BuildTutte(const MyMesh & mesh)方法为最终的构建代码，之所以这样命名是因为，上世纪60年代图特(Tutte)在解决图的嵌入问题时使用到了这个矩阵。

值得说明的是，SparseMatrix本身不支持线程安全地插入，因此在这里使用到了一个Trick，你可以查看我之前的这篇文章来了解这个Trick。

余切拉普拉斯(Cotangent Laplacian)

余切形式的拉普拉斯在三角网格上使用情形更多，在开头提出的诸多应用里，除了图特映射外，其余方法使用的拉普拉斯都为余切形式的拉普拉斯。这里的拉普拉斯不再与嵌入方式无关。它的定义为：

$L_{ij}=\left\{\begin{matrix} -\cot \alpha - \cot \beta &, e_{ij} \in E \\ \sum_{e_{ik} \in E}^{ }-L_{ik} & , i=j \\ 0 & , otherwise \end{matrix}\right.$

可以看到，除了 $e_{ij} \in E$ 的定义外，其余定义和简单拉普拉斯都是完全相同的。其中， $\cot \alpha$ 和 $\cot \beta$ 分别表示的边 $e_{ij}$ 的2个对角的余切值，如果是边界则只取1个。对应的C++代码只需要在上面的代码中添加两个函数：

//更新余切形式的拉普拉斯矩阵
template<int IsRowMajor>
void UpdateHarmonic(const MyMesh & mesh, Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor> & A)
{
	//计算非对角线元素
#pragma omp parallel for schedule(dynamic, 100)
	for (int i = 0; i < mesh.n_edges(); i++) {
		//      v1
		//    //||\\
		//   v2 || v3    
		//    \\||//
		//      v0
		const auto eh = mesh.edge_handle(i);
		const auto he01 = mesh.halfedge_handle(eh, 0);
		const auto he10 = mesh.opposite_halfedge_handle(he01);;
		const auto he12 = mesh.next_halfedge_handle(he01);
		const auto he03 = mesh.next_halfedge_handle(he10);
		const auto v0 = mesh.from_vertex_handle(he01);
		const auto v1 = mesh.to_vertex_handle(he01);
		const auto v2 = mesh.to_vertex_handle(he12);
		const auto v3 = mesh.to_vertex_handle(he03);
		const auto p0 = mesh.point(v0);
		const auto p1 = mesh.point(v1);
		const auto p2 = mesh.point(v2);
		const auto p3 = mesh.point(v3);

		MyMesh::Scalar cot = 0.0;
		if (!mesh.is_boundary(he01)) {
			const MyMesh::Point v21 = p2 - p1;
			const MyMesh::Point v20 = p2 - p0;
			cot += OpenMesh::dot(v21, v20) / OpenMesh::cross(v21, v20).norm();
		}
		if (!mesh.is_boundary(he10)) {
			const MyMesh::Point v31 = p3 - p1;
			const MyMesh::Point v30 = p3 - p0;
			cot += OpenMesh::dot(v31, v30) / OpenMesh::cross(v31, v30).norm();
		}

		MyMesh::Scalar value = -cot;
		A.coeffRef(v0.idx(), v1.idx()) = value;
		A.coeffRef(v1.idx(), v0.idx()) = value;
	}

	//计算对角线元素
	UpdateDiagonal(A);
}

//构建余切形式的拉普拉斯矩阵
template<int IsRowMajor>
Eigen::SparseMatrix<MyMesh::Scalar, IsRowMajor> BuildHarmonic(const MyMesh& mesh) {
	auto A = PrepareLaplacian<IsRowMajor>(mesh);
	UpdateHarmonic(mesh, A);
	return A;
}

上面用到了向量方法计算余切值。至于UpdateDiagonal和PrepareLaplacian方法与之前都是完全一样的，这也是为什么把这两个方法单独提取出来的原因。至于BuildHarmonic和UpdateHarmonic的分离，是因为计算参数化时通常都是迭代求解，在网格本身拓扑不变的情况下，后续迭代只需要在前面分配的矩阵上更新值即可，而不需要再次分配空间。

其它形式

当然，三角网格的拉普拉斯矩阵还有很多其它形式，不过基本形式和前面两个例子没有太多差别，主要还是某条边 $e_{ij}$ 所对应的值的具体计算方法。在此不再枚举。