【数字几何处理】 Laplacian Smooth 源码+介绍

Energy-based formulation

$u^{\ast }=argmi{n}_{u}E(u)=agrmi{n}_u\frac{1}{2}{\int }_S(f-u{)}^2+\lambda |||▽|{|}^{2}dA$
$f:noisysignal$
$\lambda :controlstherateofsmoothing$
$u$是我们最终要求解的，其中$(f-u{)}^2$项尽量减少result和noisy signal之间的差距，项$\lambda |||▽|{|}^2$则尽量控制结果更加的圆滑
注意：
1.当noisy signal为vertex’s position时，我们对顶点的position进行了平滑。除此之外我们还可以选择其他signal平滑，比如vertex’s normal。
2.energy-based formulation 和下文提到的flow-based formulation是等价的。

Flow-based formulation

$\frac{\partial u}{\partial t}=\lambda \Delta u$
$\Delta u:$the laplacian of signal u
$\lambda ：$controls the rate of smooth
显式表达：$u^{t+1}=u^t+\lambda \Delta u^t$
隐式表达：$u^{t+1}=u^t+\lambda \Delta u^{t+1}\to u^t=(I-\lambda \Delta u)u^{t+1}$

Discrete Laplacian

离散拉普拉斯算子的推导有很多种方式，由于式子太多，敲起来比较麻烦，这里省略了，直接给出结果：
$\Delta f(v_i)=\frac{1}{2A_i}{\Sigma }_{j\in {\Omega }_i}(cot{\alpha }_{ij}+cot{\beta }_{ij})(f_j-f_i)$

$\Delta f(v_i)=\frac{1}{2A_i}[-{\Sigma }_{j\in {\Omega }_i}(cot{\alpha }_{ij}+cot{\beta }_{ij})f_i+{\Sigma }_{j\in {\Omega }_i}(cot{\alpha }_{ij}+cot{\beta }_{ij})f_j]$

$\Delta f(v_i)=\left[\begin{array}{cc}-{\Sigma }_{j\in {\Omega }_i}(cot{\alpha }_{ij}+cot{\beta }_{ij})f_i& {\Sigma }_{j\in {\Omega }_i}(cot{\alpha }_{ij}+cot{\beta }_{ij})f_j\end{array}\right]\left[\begin{array}{c}{f}_i\\ f_j\end{array}\right]$

$\left[\begin{array}{c}\Delta f_i\\ \Delta f_j\\ \Delta f_k\\ ⋮\end{array}\right]=\left[\begin{array}{cccc}\frac{1}{A_i}& & & \cdots \\ & \frac{1}{A_j}& & \cdots \\ & & \frac{1}{A_k}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{cccc}-{\Sigma }_{l\in {\Omega }_i}(cot{\alpha }_{il}+cot{\beta }_{il})& cot{\alpha }_{ij}+cot{\beta }_{ij}& cot{\alpha }_{ik}+cot{\beta }_{ik}& \cdots \\ & & & \cdots \\ & & & \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{c}{f}_i\\ f_j\\ f_k\\ ⋮\end{array}\right]$
上面右边的式子中第一个矩阵为$M^{-1}$，其中M为mass matrix，第二个矩阵为cotangent matrix，最后一个矩阵为signal matrix。

求解

隐式表示的矩阵形式：
$u^t=(I-\lambda M^{-1}L)u^{t+1}$
$M{u}^t=(M-\lambda L)u^{t+1}.$
$A=M-\lambda L$
$b=M{u}^t$

源码

#include "massmatrix.h"
#include "doublearea.h"
void cgpl::massmatrix(
  const Eigen::MatrixXd & l,
  const Eigen::MatrixXi & F,
  Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
{
	Eigen::VectorXd dblA(F.rows());
	cgpl::doublearea(l, dblA);

	Eigen::VectorXd areas = Eigen::VectorXd::Zero(F.maxCoeff() + 1);

	for (int i = 0; i < F.rows(); i++) 
	{
		Eigen::Vector3i f = F.row(i);
		double area = dblA[i];

		for (int j = 0; j < f.size(); j++) 
		{
			int v = f[j];
			areas[v] += area;

		}
	}

	M.resize(F.maxCoeff() + 1);
	M.diagonal() = 1 / (double)6 * areas;
}

计算mass matrix，每个顶点的Local average region近似为周围三角形的面积除以3

#include "cotmatrix.h"
#include <iostream>
#include <vector>
#include <cmath>

typedef Eigen::Triplet<double> T;

// Given edge lengths of a triangle ABC
// Return the cotan of angle A
double my_cot(double a, double b, double c) {
	// Heron's Formula
	double s = (a + b + c) / 2.0;
	double area = sqrt(s * (s - a) * (s - b) * (s - c));

	// Law of Sine and Cosine
	double sinA = 2.0 * area / (b * c);
	double cosA = (a*a - b*b - c*c) / (-2.0 * b * c);

	return cosA / sinA;
}
void cgpl::cotmatrix(
	const Eigen::MatrixXd & l,
	const Eigen::MatrixXi & F,
	Eigen::SparseMatrix<double> & L)
{
	int nV = F.maxCoeff() + 1;
	std::vector<T> tripletList;

	for(int fi = 0; fi < F.rows(); fi++) {
		for(int e = 0; e < 3; e++) {
			int i = F(fi, e); 
			int j = F(fi, (e + 1) % 3);

			double a = l(fi, (e + 2) % 3);
			double b = l(fi, (e + 0) % 3);
			double c = l(fi, (e + 1) % 3);

			double cotA = my_cot(a, b, c);

			double Lij = 0.5 * cotA;

			tripletList.push_back(T(i, j, Lij)); 
			tripletList.push_back(T(j, i, Lij));

			// Diagonal
			tripletList.push_back(T(i, i, -Lij));
			tripletList.push_back(T(j, j, -Lij));
		}
	}

	L.resize(nV, nV);
	L.setFromTriplets(tripletList.begin(), tripletList.end());
}

计算L矩阵

#include "smooth.h"
#include "cotmatrix.h"
#include "massmatrix.h"
#include "edge_lengths.h"
void cgpl::smooth(
	const Eigen::MatrixXd & V,
	const Eigen::MatrixXi & F,
	const Eigen::MatrixXd & G,
	double lambda,
	Eigen::MatrixXd & U)
{
	Eigen::MatrixXd l;
	cgpl::edge_lengths(V, F, l);

	Eigen::SparseMatrix<double> L;
	cotmatrix(l, F, L);

	Eigen::DiagonalMatrix<double,Eigen::Dynamic> M;
	massmatrix(l, F, M);

	// M * G = (M - lambda * L) * U 
	Eigen::SparseMatrix<double> A = -lambda * L;
	for(int i = 0; i < M.rows(); i++) {
		A.coeffRef(i, i) += M.diagonal()[i];
	}

	Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver(A);
	U = solver.solve(M * G);
}

求解方程组

结果

没有smooth

smooth 3次