因为课题需要，除了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{//确定点ABCPVERT 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这个线性方程组，求解结果就是每个顶点对应的新坐标值，最后将其赋值给变形点就完成了整个变形。
变形效果如图，变形前：
变形后：
效果还是不错的。
