因为课题需要，除了RBF还做了一个Laplace网格变形，其他大佬已经把原理写的很详细了，我就简单介绍一下公式，主要还是写写实现过程。过程同样参考了大佬的部分代码，而且实现的时候刚开始敲代码不久，所以有点乱QAQ。
首先，计算离散拉普拉斯坐标，网格上的点vi的拉普拉斯坐标δi为：
\[\delta_{i}=v_{i}-\sum_{j \in N(i)} \omega_{i j} v_{j} \]式中，ωij为点vj相对vi的权值，如图所示，点vj就是点vi的邻接点。本文实现的权值有均匀权值和余切权值。
均匀权值的计算公式为：
\[\mathit{\omega {\tiny ij} }=\frac{1}{d{\tiny i} } \]di就是邻接顶点的数量，代码也相对于余切权值简单不少，但是变形效果较差。
余切权值的计算公式为：
\[\mathit{\omega {\tiny ij} }=\frac{\cot \alpha {\tiny ij} +\cot \beta {\tiny ij} }{2 } \]αij、βij其实就是边ij的两个对角，如图所示：
然后，根据输入的待变形顶点点集，移动点点集和固定点点集，构建拉普拉斯矩阵，然后通过求解一个线性方程组就可以获取变形后的顶点位置：
\[Lv_{}^{\prime } = \Delta \]其中，L就是拉普拉斯矩阵，v为变形后的顶点，Δ为拉普拉斯坐标矩阵，剩下的流程在代码里说明。
首先是构建权值矩阵：
void Laplace::calculateCotValue(Vector3d A, Vector3d B, Vector3d C,double &cotValue)
{
Eigen::Vector3d Vab = B - A;
Eigen::Vector3d Vac = C - A;
double cosValue = Vab.dot(Vac) / (Vab.norm() * Vac.norm());
double sinValue = (Vab.cross(Vac)).norm() / (Vab.norm() * Vac.norm());
cotValue = cosValue / sinValue;
}
void Laplace::baseCotWeight(STLModel* p_STL)
{
double smoothingFactor = 0.0001; // 光顺因子
/*-------------------------*/
//建立邻接表
int Vert_num = this->m_DeforPoint.size();//参与变形的数量
int Move_num = this->m_MoveIndex.size();//移动锚点的数量
int Fix_num = this->m_FixPoint.size();//固定锚点的数量
unordered_set<BaseStruct::PVERT> Jud;//用于判断变形点是否为所选择的点
for (auto it : this->m_DeforPoint)
{
Jud.insert(it);
}
Weight.resize(Vert_num, Vert_num);
L.resize(Vert_num, Vert_num);
Ls.resize(Vert_num, Vert_num);
/*------------------邻接点赋值为wij，i点处赋值为1，其他非邻接赋值为0------------------*/
typedef Eigen::Triplet<double> T_Cot;//创造三元组记录余切值
std::vector<T_Cot> tripletList;
for (int i = 0;i < Vert_num;i++)
{
BaseStruct::PVERT temp_vert = this->m_DeforPoint[i];
BaseStruct::PHEDGE start = temp_vert->pHEdgeOut->pHEdgeNext;
BaseStruct::PHEDGE now = start;
double sum_weight = 0.0; //权值总和
map<unsigned int, double>vWeight;
do
{
//确定点ABC
PVERT Target = now->pHEdgePair->pVertEnd;
double cotValue;double cotValue1;double cotValue2;
Eigen::Vector3d A(now->pVertEnd->x, now->pVertEnd->y, now->pVertEnd->z);
Eigen::Vector3d B(now->pHEdgePair->pVertEnd->x, now->pHEdgePair->pVertEnd->y, now->pHEdgePair->pVertEnd->z);
Eigen::Vector3d C(now->pHEdgeNext->pVertEnd->x, now->pHEdgeNext->pVertEnd->y, now->pHEdgeNext->pVertEnd->z);
Eigen::Vector3d D(now->pHEdgePair->pHEdgeNext->pVertEnd->x, now->pHEdgePair->pHEdgeNext->pVertEnd->y, now->pHEdgePair->pHEdgeNext->pVertEnd->z);
calculateCotValue(B, A, C, cotValue1);
calculateCotValue(B, A, D, cotValue2);
cotValue = (0.5) * (cotValue1 + cotValue2 + epsilon);
if (Jud.find(Target) != Jud.end())
{
vWeight.insert(pair<unsigned int, double>(now->pHEdgePair->pVertEnd->index, cotValue));
sum_weight += cotValue;
}
now = now->pHEdgePair->pHEdgePre;
} while (start != now);
for (map< unsigned int, double>::iterator it = vWeight.begin();it!=vWeight.end();it++)
{
double cotvalue = it->second / sum_weight;
tripletList.push_back(T_Cot(temp_vert->index, it->first, -cotvalue));
}
tripletList.push_back(T_Cot(temp_vert->index, temp_vert->index, 1));
}
Weight.setFromTriplets(tripletList.begin(), tripletList.end());
Ls = L = Weight;
Ls.conservativeResize(Vert_num + Fix_num + Move_num, Vert_num);//添加锚点
/*-------------------------------------------------------*/
//计算Ls和LsT、LsTLs
/*-------------------------------------------------------*/
for (int a = 0;a < this->m_FixPoint.size();a++)
{
Ls.coeffRef(a + Vert_num, this->m_FixPoint[a]->index) = 1;
}
int a = 0;
for (auto it : m_MoveIndex)
{
Ls.coeffRef(a +Vert_num+Fix_num, it) = 1;
a++;
}
LsT = Ls.transpose();
LsTLs = LsT * Ls;
这段代码中，第二个for循环就是计算拉普拉斯矩阵的主要部分，这部分我用到的方法基于半边结构实现的，而且写的很乱，但是逻辑就是前面介绍的余切权值计算方法。计算完网格的拉普拉斯矩阵L后，再构建加入固定点和变形点的拉普拉斯矩阵Ls，最后求转置再相乘的原因是为了让其可逆，以便等等能够求解变形后的顶点坐标。
void Laplace::BuildLsTb(STLModel* p_STL)
{
/*----------顶点、移动点、固定点数量-------------*/
const int move_anchors_num = Move_anchor_coor.size();
const int Vext_num = this->m_DeforPoint.size();
const int Fix_num = this->m_FixPoint.size();//固定锚点的数量
/*-----------计算拉普拉斯坐标B-------------------*/
vx.resize(Vext_num);vy.resize(Vext_num);vz.resize(Vext_num);//每一个顶点x,y,z坐标
for (int i = 0;i < Vext_num;i++)
{
vx.coeffRef(i) = this->m_DeforPoint[i]->x;
vy.coeffRef(i) = this->m_DeforPoint[i]->y;
vz.coeffRef(i) = this->m_DeforPoint[i]->z;
}
bx = L * vx;
by = L * vy;
bz = L * vz;//变形前的拉普拉斯坐标值
bx.conservativeResize(Vext_num + Fix_num + move_anchors_num);
by.conservativeResize(Vext_num + Fix_num + move_anchors_num);
bz.conservativeResize(Vext_num + Fix_num + move_anchors_num);
// 用形变后坐标对移动锚点坐标进行赋值
int i = 0;
for (auto it : this->m_FixPoint)
{
bx.coeffRef(Vext_num + i) = it->x;
by.coeffRef(Vext_num + i) = it->y;
bz.coeffRef(Vext_num + i) = it->z;
i++;
}
for (int i = 0;i < move_anchors_num;i++)
{
bx.coeffRef(Vext_num + Fix_num + i) = Move_anchor_coor[i].x;
by.coeffRef(Vext_num + Fix_num + i) = Move_anchor_coor[i].y;
bz.coeffRef(Vext_num + Fix_num + i) = Move_anchor_coor[i].z;
}
LsTbx = LsT * bx;
LsTby = LsT * by;
LsTbz = LsT * bz;
}
这一段代码目的很明确，就是计算网格顶点变形前的拉普拉斯坐标，并构建坐标矩阵。
void Laplace :: LaplaceDeformation(STLModel* p_STL)
{
LsTLs.makeCompressed();
SparseLU<SparseMatrix<double>> LU;
LU.analyzePattern(LsTLs);
LU.factorize(LsTLs);
LU.info();
vx_New = LU.solve(LsTbx);
vy_New = LU.solve(LsTby);
vz_New = LU.solve(LsTbz);
}
然后计算变形后的坐标值，这里我把x，y，z三个值分别拿出来进行计算，根据计算公式，最后求解Ax=b这个线性方程组，求解结果就是每个顶点对应的新坐标值，最后将其赋值给变形点就完成了整个变形。
变形效果如图，变形前：
变形后：
效果还是不错的。
标签：变形,拉普拉斯,网格,pVertEnd,int,num,Vext,now From： https://www.cnblogs.com/nanjofish/p/18277158