【数字几何处理】Deformation：Laplacian-based energy&&As-rigid-as-possible 源码+介绍

Laplacian-based energy

原理：
$mi{n}_x{\int }_S||\Delta x-\Delta \hat{x}|{|}^{2}dA$

$d=x-\hat{x}$

=>$mi{n}_d{\int }_S||\Delta d|{|}^{2}dA$

=>$mi{n}_D{D}^T{L}^T{M}^{-T}M{M}^{-1}LD$

$M:MassMatrix$
$L:cotMatrix$

=>$mi{n}_D{D}^T{L}^T{M}^{-1}LD$

=>$mi{n}_D{D}^{T}QD$

$D_{handles}$已知

=>$mi{n}_D{D}^{T}QD$ subject to $D_{handles}=D_{kownValue}$

每次拖动鼠标时，变化的只是$D_{handles}$，所以可以先计算出Q，在实时计算时实时更新$D_{handles}$即可，算出D后，$V_{new}=V_{old}+D$

#include "biharmonic_precompute.h"

#include "cotmatrix.h"
#include "massmatrix.h"
void biharmonic_precompute(
	const Eigen::MatrixXd & V,
	const Eigen::MatrixXi & F,
	const Eigen::VectorXi & b,
	Eigen::SparseMatrix<double> & QQ)
{
	Eigen::SparseMatrix<double> L;
	cgpl::cotmatrix(V, F, L);

	Eigen::SparseMatrix<double> M;
	cgpl::massmatrix(V, F, M);

	Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver(M);
	Eigen::SparseMatrix<double> M_inv_L = solver.solve(L).sparseView();
	Eigen::SparseMatrix<double> Q = L.transpose() * M_inv_L; // L^T M^-1 L

	typedef Eigen::Triplet<double>T;
	std::vector<T>tripleList;

	std::vector<bool> Known(V.rows(), false);
	for (int i = 0; i < b.rows(); i++)
	{
		Known[b(i)] = true;
		tripleList.push_back(T(b(i), b(i), 1));
	}


	for (int i = 0; i < Q.outerSize(); i++)
	{
		for (Eigen::SparseMatrix<double>::InnerIterator it(Q, i); it; ++it)
		{
			int row = it.row();
			int col = it.col();
			if (Known[row] == false)
			{
				tripleList.push_back(T(row, col, it.value()));
			}
		}
	}

	QQ.resize(V.rows(), V.rows());
	QQ.setFromTriplets(tripleList.begin(), tripleList.end());
}

#include "biharmonic_solve.h"
#include<Eigen/SparseLU>
void biharmonic_solve(
	const Eigen::SparseMatrix<double> & QQ,
	const Eigen::VectorXi & b,
	const Eigen::MatrixXd & bc,
	Eigen::MatrixXd & D)
{
	D.resize(QQ.cols(), 3);
	Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>solver;

	Eigen::VectorXd dx = Eigen::VectorXd::Zero(QQ.cols());
	Eigen::VectorXd dy = Eigen::VectorXd::Zero(QQ.cols());
	Eigen::VectorXd dz = Eigen::VectorXd::Zero(QQ.cols());
	for (int i = 0; i < b.rows(); i++)
	{
		dx(b(i)) = bc(i, 0);
		dy(b(i)) = bc(i, 1);
		dz(b(i)) = bc(i, 2);
	}

	solver.compute(QQ);
	Eigen::VectorXd Dx=solver.solve(dx);
	Eigen::VectorXd Dy=solver.solve(dy);
	Eigen::VectorXd Dz=solver.solve(dz);

	D << Dx, Dy, Dz;
}

As-rigid-as-possible

原理：
Rigid Transformation: $x_i=R\widehat{x_i}+t$
$\frac{1}{2}{\int }_{\Omega }\Vert ▽x-▽\hat{x}{\Vert }^{2}dA$

=>$\frac{1}{2}{\Sigma }_{f\in F}{\Sigma }_{ij\in f}{c}_{ij}\Vert (v_i-v_j)-({\hat{v}}_i-{\hat{v}}_j)\Vert ²$

=>$\frac{1}{6}{\Sigma }_{k=1}^N{\Sigma }_{ij\in F(K)}{c}_{ij}\Vert (v_i-v_j)-({\hat{v}}_i-{\hat{v}}_j)\Vert ²$

=>$\frac{1}{6}{\Sigma }_{k=1}^N{\Sigma }_{ij\in F(K)}{c}_{ij}(v_i-v_j{)}^T(v_i-v_j)+\frac{1}{6}{\Sigma }_{k=1}^N{\Sigma }_{ij\in F(K)}{c}_{ij}(v_i-v_j{)}^T{R}_k({\hat{v}}_i-{\hat{v}}_j)$

=>$tr(V^{T}LV)+tr(V^{T}KR)$

Local step（固定V优化R）：

$mi{n}_{R}tr(V^{T}KR)$

$C^T=V^{T}K$

用到的一个结论：

$C^T=U\Sigma V^T$

$R=U{V}^T$时$tr(V^{T}KR)$有最优解

Global step（固定R优化V）

$mi{n}_{V}tr(V^{T}LV)+tr(V^{T}B)$

#include "arap_precompute.h"
#include "cotmatrix.h"
typedef Eigen::Triplet<double> T;
void arap_precompute(
	const Eigen::MatrixXd & V,
	const Eigen::MatrixXi & F,
	const Eigen::VectorXi & b,
	Eigen::SparseMatrix<double> & Q,
	Eigen::SparseMatrix<double> & K)
{
	Eigen::SparseMatrix<double> L;
	cgpl::cotmatrix(V, F, L);

	std::vector<T>tripleList;

	std::vector<bool> Known(V.rows(), false);
	for (int i = 0; i < b.rows(); i++)
	{
		Known[b(i)] = true;
		tripleList.push_back(T(b(i), b(i), 1));
	}


	for (int i = 0; i < L.outerSize(); i++)
	{
		for (Eigen::SparseMatrix<double>::InnerIterator it(L, i); it; ++it)
		{
			int row = it.row();
			int col = it.col();
			if (Known[row] == false)
			{
				tripleList.push_back(T(row, col, it.value()));
			}
		}
	}

	Q.resize(V.rows(), V.rows());
	Q.setFromTriplets(tripleList.begin(), tripleList.end());

	int nV = V.rows();
	std::vector<T> tripletList;

	for(int f = 0; f < F.rows(); f++) {
		for(int a = 0; a < 3; a++) {
			int i = F(f, (a + 1) % 3); 
			int j = F(f, (a + 2) % 3);

			Eigen::Vector3d diff = V.row(i) - V.row(j);
			Eigen::Vector3d eij = L.coeff(i, j) * diff / 6.0;
			for(int ki = 0; ki < 3; ki++) {
				int k = F(f, (a + ki) % 3);
				for(int B = 0; B < 3; B++) {
					tripletList.push_back(T(i, 3 * k + B, eij[B])); 
					tripletList.push_back(T(j, 3 * k + B, -eij[B]));
				}
			}
		}
	}

	K.resize(nV, 3 * nV);
	K.setFromTriplets(tripletList.begin(), tripletList.end());
}

#include "arap_single_iteration.h"
void polar_svd3x3(const Eigen::Matrix3d &C,Eigen::Matrix3d &R)
{
	Eigen::JacobiSVD<Eigen::MatrixXd> svd(C, Eigen::ComputeThinU | Eigen::ComputeThinV);
	Eigen::Matrix3d V = svd.matrixV();
	Eigen::Matrix3d U = svd.matrixU();
	R = U* V.transpose() ;
}
void arap_single_iteration(
	const Eigen::SparseMatrix<double> & Q,
	const Eigen::VectorXi & b,
	const Eigen::SparseMatrix<double> & K,
	const Eigen::MatrixXd & bc,
	Eigen::MatrixXd & U)
{
	Eigen::MatrixXd C = K.transpose() * U;
	Eigen::MatrixXd R(C.rows(), C.cols());

	for(int k = 0; k < Q.cols(); k++) {
		Eigen::Matrix3d C_k = C.block(k * 3, 0, 3, 3);
		Eigen::Matrix3d R_k;
		C_k.normalize();
		polar_svd3x3(C_k, R_k);
		R.block(k * 3, 0, 3, 3) = R_k;
	}

	Eigen::MatrixXd B= K * R;
	B = -B;

	U.resize(Q.cols(), 3);
	Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>solver;

	for (int i = 0; i < b.rows(); i++)
	{
		B(b(i),0) = bc(i, 0);
		B(b(i),1) = bc(i, 1);
		B(b(i),2) = bc(i, 2);
	}

	solver.compute(Q);
	Eigen::VectorXd X=solver.solve(B.col(0));
	Eigen::VectorXd Y=solver.solve(B.col(1));
	Eigen::VectorXd Z=solver.solve(B.col(2));

	U << X, Y, Z;
}

结果

Laplacian-based energy：

As-rigid-as-possible