1.设置网格顶点局部标架
定义顶点 Vi 的局部标架 Fi = (ei1, ei2, ni)，如图
三维空间中的任意向量 A 可用局部标架表示为 A = λ1e1 + λ2e2 + λ3n;
2.求取矩阵 T
两两局部标架之间有旋转矩阵 T（因为是局部标架而非局部坐标系，所以没有平移，不适用合同变换），如下
每个顶点的标架矩阵 F 都是3 x 3的矩阵，设Fi = (ei1, ei2, ni)和 Fj = (ej1, ej2, nj)为顶点Vi 和 Vj的标架矩阵，Vj 为Vi 的一环邻域顶点，如图
则对于旋转矩阵 Tij，有 Fj - TijFi = 0; 调整一下为 Tij 等于 Fj 乘以 Fi的逆。当世界坐标系发生变换时，标架 Fi 和 Fj 同样发生变化，但是矩阵 Tij 时不会变化的。
此为网格的内在属性。
3.求取更新后的局部标架
和之前那篇，求取拉普拉斯矩阵差不多。只是网格顶点变化为 3 x 3 的矩阵而已，如下
关于矩阵 G ，其中每一行的分量是为顶点 Vi 对于其一阶邻域顶点 Vj 的旋转矩阵 Tji ，每一列的分量是为 一阶邻域顶点 Vj 到 顶点 Vi 的旋转矩阵 Tji
关于矩阵 H ，和拉普拉斯矩阵中的约束调参顶点一样，但是要修改为 3 x 3的单位矩阵 i。 当 R矩阵中的约束顶点标架发生变化时，便可以求取新的顶点标架 F ‘ 。
4. 求取旋转不变网格之后的拉普拉斯网格形变
和之前求取拉普拉斯坐标δ一样，只是需要将全局拉普拉斯坐标转换为局部标架下的局部拉普拉斯坐标δi
, 
当步骤3求取到每个顶点的新标架F之后，求取新的局部拉普拉斯坐标，进而求得拉普拉斯坐标。之后就可以和求取拉普拉斯算子矩阵方程一样，解稀疏线性方程组。
构建局部标架以及原始旋转矩阵T和局部拉普拉斯坐标
1 void mainNode::InitAxesAndMatrixSelf() 2 { 3 //先构建局部坐标系 4 for (int i = 0; i < m_vecAllVertex.size(); i++) 5 { 6 pVERTEX pVertex = m_vecAllVertex.at(i); 7 std::vector<pTRIANGLE> vectorTri = pVertex->vecNeighborTri; 8 osg::Vec3 vertexNormal, E1, E2; 9 vertexNormal = E1 = E2 = osg::Vec3(0.0, 0.0, 0.0); 10 int sizeOfTri = vectorTri.size(); 11 for (int j = 0; j < sizeOfTri; j++) 12 { 13 pTRIANGLE pTri = vectorTri.at(j); 14 pVERTEX pA = pTri->pA; 15 pVERTEX pB = pTri->pB; 16 pVERTEX pC = pTri->pC; 17 osg::Vec3 BA = pA->pos - pB->pos; 18 osg::Vec3 BC = pC->pos - pB->pos; 19 osg::Vec3 normal = BA ^ BC; 20 normal.normalize(); 21 vertexNormal += normal; 22 } 23 vertexNormal /= sizeOfTri; 24 E1 = osg::Vec3( - vertexNormal.y(), vertexNormal.x(), 0.0); 25 E2 = vertexNormal ^ E1; 26 vertexNormal.normalize(); 27 E1.normalize(); 28 E2.normalize(); 29 pVertex->E1 = E1; 30 pVertex->E2 = E2; 31 pVertex->N = vertexNormal; 32 pVertex->axesSelf.resize(3, 3); 33 std::vector<Eigen::Triplet<float>> vectorTriplet; 34 vectorTriplet.push_back(Eigen::Triplet<float>(0, 0, E1.x())); 35 vectorTriplet.push_back(Eigen::Triplet<float>(1, 0, E1.y())); 36 vectorTriplet.push_back(Eigen::Triplet<float>(2, 0, E1.z())); 37 38 vectorTriplet.push_back(Eigen::Triplet<float>(0, 1, E2.x())); 39 vectorTriplet.push_back(Eigen::Triplet<float>(1, 1, E2.y())); 40 vectorTriplet.push_back(Eigen::Triplet<float>(2, 1, E2.z())); 41 42 vectorTriplet.push_back(Eigen::Triplet<float>(0, 2, vertexNormal.x())); 43 vectorTriplet.push_back(Eigen::Triplet<float>(1, 2, vertexNormal.y())); 44 vectorTriplet.push_back(Eigen::Triplet<float>(2, 2, vertexNormal.z())); 45 46 Eigen::SparseMatrix<float> matrixTemp(3, 3); 47 matrixTemp.setFromTriplets(vectorTriplet.begin(), vectorTriplet.end()); 48 pVertex->axesSelf = matrixTemp; 49 //构建相对局部坐标系的拉普拉斯坐标长度分量 50 pVertex->relativeX = pVertex->lplsPos * E1; 51 pVertex->relativeY = pVertex->lplsPos * E2; 52 pVertex->relativeZ = pVertex->lplsPos * vertexNormal; 53 } 54 //构建周围一阶顶点到此顶点的转换矩阵 55 for (int i = 0; i < m_vecAllVertex.size(); i++) 56 { 57 pVERTEX pVertex = m_vecAllVertex.at(i); 58 Eigen::Matrix3f axesSelf = pVertex->axesSelf; 59 //Eigen::Matrix3f matrixSelf = axesSelf; 60 std::vector<pVERTEX> vecNeighborVertex = pVertex->vecNeighborVertex; 61 for (int j = 0; j < vecNeighborVertex.size(); j++) 62 { 63 pVERTEX pNeighbor = vecNeighborVertex.at(j); 64 Eigen::Matrix3f axesOther = pNeighbor->axesSelf; 65 //std::cout << axesOther << std::endl; 66 Eigen::Matrix3f axesOtherInverse = axesOther.inverse(); 67 //std::cout << axesOtherInverse << std::endl; 68 // Ax = b; T * axesOther = axesSelf; 69 Eigen::Matrix3f matrixT = axesSelf * axesOtherInverse; 70 pVertex->mapMatrixOtherToSelf.insert(std::pair<int, Eigen::Matrix3f>(pNeighbor->id, matrixT)); 71 //test 72 Eigen::Matrix3f matrixTemp = matrixT * axesOther; 73 //std::cout << matrixTemp << std::endl; 74 //std::cout << axesSelf << std::endl; 75 } 76 } 77 }
构建旋转矩阵方程组
1 void mainNode::InitNewAxesMatrix() 2 { 3 int sizeOfVertex = m_vecAllVertex.size(); 4 std::vector<Eigen::Triplet<float> > vectorTriplet; 5 for (int i = 0; i < sizeOfVertex; i++) 6 { 7 pVERTEX pVertex = m_vecAllVertex.at(i); 8 int vertexID = pVertex->id; 9 float sizeOfNeighborVertex = pVertex->vecNeighborVertex.size(); 10 11 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3, vertexID * 3, - sizeOfNeighborVertex)); 12 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 1, vertexID * 3 + 1, - sizeOfNeighborVertex)); 13 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 2, vertexID * 3 + 2, - sizeOfNeighborVertex)); 14 15 std::map<int, Eigen::Matrix3f>::iterator itorOfT; 16 for(itorOfT = pVertex->mapMatrixOtherToSelf.begin(); itorOfT != pVertex->mapMatrixOtherToSelf.end(); itorOfT++) 17 { 18 int vertexNeighborID = itorOfT->first; 19 Eigen::Matrix3f T = itorOfT->second; 20 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3, vertexNeighborID * 3, T(0, 0))); 21 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3, vertexNeighborID * 3 + 1, T(0, 1))); 22 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3, vertexNeighborID * 3 + 2, T(0, 2))); 23 24 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 1, vertexNeighborID * 3, T(1, 0))); 25 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 1, vertexNeighborID * 3 + 1, T(1, 1))); 26 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 1, vertexNeighborID * 3 + 2, T(1, 2))); 27 28 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 2, vertexNeighborID * 3, T(2, 0))); 29 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 2, vertexNeighborID * 3 + 1, T(2, 1))); 30 vectorTriplet.push_back(Eigen::Triplet<float>(vertexID * 3 + 2, vertexNeighborID * 3 + 2, T(2, 2))); 31 } 32 } 33 //添加约束顶点;头一圈和尾一圈,可自己调节 34 for (int i = 0; i < m_iAroundNumber * 2; i++) 35 { 36 pVERTEX pVertex; 37 int vertexID; 38 float value = 1.0; 39 if (i < m_iAroundNumber) 40 { 41 pVertex = m_vecAllVertex.at(i); 42 vertexID = pVertex->id; 43 } 44 else 45 { 46 pVertex = m_vecAllVertex.at(i + m_iAroundNumber * (m_iHeightNumber - 2)); 47 vertexID = pVertex->id; 48 } 49 vectorTriplet.push_back(Eigen::Triplet<float>(sizeOfVertex * 3 + i * 3, vertexID * 3, value)); 50 vectorTriplet.push_back(Eigen::Triplet<float>(sizeOfVertex * 3 + i * 3 + 1, vertexID * 3 + 1, value)); 51 vectorTriplet.push_back(Eigen::Triplet<float>(sizeOfVertex * 3 + i * 3 + 2, vertexID * 3 + 2, value)); 52 } 53 54 m_NewAxesMatrix.resize(sizeOfVertex * 3 + (m_iAroundNumber * 2 * 3), sizeOfVertex * 3); 55 m_NewAxesMatrix.setFromTriplets(vectorTriplet.begin(), vectorTriplet.end()); 56 }
求取新的局部标架并获得新的拉普拉斯坐标
1 void mainNode::CalculateNewAxes() 2 { 3 Eigen::SparseMatrix<float> matrixNewAxesTrans = m_NewAxesMatrix.transpose(); 4 Eigen::SparseMatrix<float> matrixRes = matrixNewAxesTrans * m_NewAxesMatrix; 5 Eigen::SimplicialCholesky<Eigen::SparseMatrix<float> > matrixCholesky(matrixRes); 6 7 std::vector<Eigen::Triplet<float> > vectorTriplet; 8 int sizeOfVertex = m_vecAllVertex.size(); 9 //方程式等号右边 10 Eigen::VectorXf vectorX, vectorY, vectorZ; 11 vectorX.resize(sizeOfVertex * 3 + m_iAroundNumber * 2 * 3); 12 vectorX.setZero(); 13 vectorY.resize(sizeOfVertex * 3 + m_iAroundNumber * 2 * 3); 14 vectorY.setZero(); 15 vectorZ.resize(sizeOfVertex * 3 + m_iAroundNumber * 2 * 3); 16 vectorZ.setZero(); 17 for (int i = 0; i < m_iAroundNumber * 2; i++) 18 { 19 pVERTEX pVertex; 20 int vertexID; 21 Eigen::Matrix3f matrixSelf; 22 float value = 1.0; 23 if (i < m_iAroundNumber) 24 { 25 pVertex = m_vecAllVertex.at(i); 26 } 27 else 28 { 29 pVertex = m_vecAllVertex.at(i + m_iAroundNumber * (m_iHeightNumber - 2)); 30 } 31 vertexID = pVertex->id; 32 matrixSelf = pVertex->axesSelf; 33 34 vectorX[sizeOfVertex * 3 + i * 3] = matrixSelf(0, 0); 35 vectorX[sizeOfVertex * 3 + i * 3 + 1] = matrixSelf(1, 0); 36 vectorX[sizeOfVertex * 3 + i * 3 + 2] = matrixSelf(2, 0); 37 38 vectorY[sizeOfVertex * 3 + i * 3] = matrixSelf(0, 1); 39 vectorY[sizeOfVertex * 3 + i * 3 + 1] = matrixSelf(1, 1); 40 vectorY[sizeOfVertex * 3 + i * 3 + 2] = matrixSelf(2, 1); 41 42 vectorZ[sizeOfVertex * 3 + i * 3] = matrixSelf(0, 2); 43 vectorZ[sizeOfVertex * 3 + i * 3 + 1] = matrixSelf(1, 2); 44 vectorZ[sizeOfVertex * 3 + i * 3 + 2] = matrixSelf(2, 2); 45 } 46 47 vectorX = matrixNewAxesTrans * vectorX; 48 vectorY = matrixNewAxesTrans * vectorY; 49 vectorZ = matrixNewAxesTrans * vectorZ; 50 51 Eigen::VectorXf vectorResX, vectorResY, vectorResZ; 52 vectorResX = matrixCholesky.solve(vectorX); 53 vectorResY = matrixCholesky.solve(vectorY); 54 vectorResZ = matrixCholesky.solve(vectorZ); 55 //得到新的标架 56 for (int i = 0; i < m_vecAllVertex.size(); i++) 57 { 58 Eigen::Matrix3f matrixSelf; 59 pVERTEX pVertex = m_vecAllVertex.at(i); 60 float x1, x2, x3, y1, y2, y3, z1, z2, z3; 61 x1 = vectorResX[i * 3]; 62 x2 = vectorResX[i * 3 + 1]; 63 x3 = vectorResX[i * 3 + 2]; 64 65 y1 = vectorResY[i * 3]; 66 y2 = vectorResY[i * 3 + 1]; 67 y3 = vectorResY[i * 3 + 2]; 68 69 z1 = vectorResZ[i * 3]; 70 z2 = vectorResZ[i * 3 + 1]; 71 z3 = vectorResZ[i * 3 + 2]; 72 73 pVertex->E1 = osg::Vec3(x1, x2, x3); 74 pVertex->E2 = osg::Vec3(y1, y2, y3); 75 pVertex->N = osg::Vec3(z1, z2, z3); 76 77 //新的拉普拉斯坐标 78 pVertex->lplsPos = (pVertex->E1 * pVertex->relativeX) + (pVertex->E2 * pVertex->relativeY) + (pVertex->N * pVertex->relativeZ); 79 } 80 }
构建拉普拉斯矩阵
1 void mainNode::InitLpls() 2 { 3 int vertexNumber = m_vecAllVertex.size(); 4 //计算顶点一阶临接顶点的权值,方式一 5 for (int i = 0; i < vertexNumber; i++) 6 { 7 pVERTEX pVertex = m_vecAllVertex.at(i); 8 int neighborNumber = pVertex->vecNeighborVertex.size(); 9 for (int j = 0; j < neighborNumber; j++) 10 { 11 pVERTEX pNeighbor = pVertex->vecNeighborVertex.at(j); 12 float weight = float(1) / float(neighborNumber); 13 pVertex->mapWeightForOther.insert(std::pair<int, float>(pNeighbor->id, weight) ); 14 pVertex->fTotalWeight = 1.0; 15 } 16 } 17 18 /* 19 //构建拉普拉斯矩阵权值，方式二 20 for (int i = 0; i < vertexNumber; i++) 21 { 22 pVERTEX pVertex = m_vecAllVertex.at(i); 23 float totalWeight = 0.0; 24 for (int j = 0; j < pVertex->vecNeighborEdge.size(); j++) 25 { 26 pEDGE pEdge = pVertex->vecNeighborEdge.at(j); 27 pVERTEX pA = pEdge->pA; 28 pVERTEX pB = pEdge->pB; 29 30 pVERTEX pTarget; 31 if (pA->id == pVertex->id) 32 pTarget = pB; 33 else 34 pTarget = pA; 35 36 std::vector<pTRIANGLE> vecTri = pEdge->vecNeighborTri; 37 pVERTEX pC = NULL; 38 float weight = 0.0; 39 for (int k = 0; k < vecTri.size(); k++) 40 { 41 pTRIANGLE pTri = vecTri.at(k); 42 pVERTEX p1 = pTri->pA; 43 pVERTEX p2 = pTri->pB; 44 pVERTEX p3 = pTri->pC; 45 if ((pA->id == p1->id && pB->id == p2->id) || (pA->id == p2->id && pB->id == p1->id)) 46 { 47 pC = p3; 48 } 49 else if ((pA->id == p1->id && pB->id == p3->id) || (pA->id == p3->id && pB->id == p1->id)) 50 { 51 pC = p2; 52 } 53 else if ((pA->id == p2->id && pB->id == p3->id) || (pA->id == p3->id && pB->id == p2->id)) 54 { 55 pC = p1; 56 } 57 //开始求取权值 58 float cotAngle = 0.0; 59 GetCotAngle(pA->pos, pB->pos, pC->pos, cotAngle); 60 weight += cotAngle; 61 } 62 weight /= 2.0; 63 pVertex->mapWeightForOther.insert(std::pair<int, float>(pTarget->id, weight) ); 64 totalWeight += weight; 65 } 66 pVertex->fTotalWeight = totalWeight; 67 } 68 */ 69 //计算拉普拉斯坐标 70 for (int i = 0; i < vertexNumber; i++) 71 { 72 pVERTEX pVertex = m_vecAllVertex.at(i); 73 osg::Vec3 pos = pVertex->pos; 74 pos = pos * pVertex->fTotalWeight; 75 osg::Vec3 otherPos = osg::Vec3(0.0, 0.0, 0.0); 76 for (int j = 0; j < pVertex->vecNeighborVertex.size(); j++) 77 { 78 pVERTEX pNeihbor = pVertex->vecNeighborVertex.at(j); 79 std::map<int, float>::iterator itor = pVertex->mapWeightForOther.find(pNeihbor->id); 80 float weight = itor->second; 81 otherPos += pNeihbor->pos * weight; 82 } 83 pVertex->lplsPos = pos - otherPos; 84 } 85 86 int count = 0; 87 std::vector<int> beginNumber(vertexNumber); 88 for (int i = 0; i < vertexNumber; i++) 89 { 90 beginNumber[i] = count; 91 count += m_vecAllVertex[i]->vecNeighborVertex.size() + 1; 92 } 93 94 std::vector<Eigen::Triplet<float> > vectorTriplet(count); 95 for (int i = 0; i < vertexNumber; i++) 96 { 97 pVERTEX pVertex = m_vecAllVertex.at(i); 98 //原始拉普拉斯矩阵 99 vectorTriplet[beginNumber[i]] = Eigen::Triplet<float>(i, i, pVertex->fTotalWeight); 100 int j = 0; 101 std::map<int, float>::iterator itor; 102 for(itor = pVertex->mapWeightForOther.begin(); itor != pVertex->mapWeightForOther.end(); itor++) 103 { 104 int neighborID = itor->first; 105 float weight = itor->second; 106 vectorTriplet[beginNumber[i] + j + 1] = Eigen::Triplet<float>(i, neighborID, -weight); 107 j++; 108 } 109 } 110 111 //头一圈和顶一圈 112 for (int i = 0; i < m_iAroundNumber * 2; i++) 113 { 114 float weight = 1.0; 115 pVERTEX pVertex; 116 if (i < m_iAroundNumber) 117 { 118 pVertex = m_vecAllVertex.at(i); 119 } 120 else 121 { 122 pVertex = m_vecAllVertex.at(i + m_iAroundNumber * 14); 123 } 124 125 int row = i + vertexNumber; 126 vectorTriplet.push_back(Eigen::Triplet<float>(row, pVertex->id, weight)); 127 } 128 129 m_Matrix.resize(vertexNumber + m_iAroundNumber * 2, vertexNumber); 130 //m_Matrix.resize(vertexNumber, vertexNumber); 131 m_Matrix.setFromTriplets(vectorTriplet.begin(), vectorTriplet.end()); 132 //std::cout << m_Matrix << std::endl; 133 }
根据新的拉普拉斯坐标求取新的顶点坐标
1 void mainNode::CalculateNewMatrixLpls() 2 { 3 Eigen::SparseMatrix<float> Matrix = m_Matrix; 4 5 Eigen::SparseMatrix<float> MatrixTranspose = Matrix.transpose(); 6 Eigen::SparseMatrix<float> MatrixLpls = MatrixTranspose * Matrix; 7 8 Eigen::VectorXf targetX, targetY, targetZ; 9 int vertexNumber = m_vecAllVertex.size(); 10 11 targetX.resize(vertexNumber + m_iAroundNumber * 2); 12 targetY.resize(vertexNumber + m_iAroundNumber * 2); 13 targetZ.resize(vertexNumber + m_iAroundNumber * 2); 14 15 //方程式等号右边 16 for (int i = 0; i < vertexNumber + m_iAroundNumber * 2; i++) 17 { 18 if (i < vertexNumber) 19 { 20 targetX[i] = m_vecAllVertex.at(i)->lplsPos[0]; 21 targetY[i] = m_vecAllVertex.at(i)->lplsPos[1]; 22 targetZ[i] = m_vecAllVertex.at(i)->lplsPos[2]; 23 } 24 else if (i < vertexNumber + m_iAroundNumber) 25 { 26 //第一层的顶点坐标分量 27 targetX[i] = m_vecAllVertex.at(i - vertexNumber)->pos.x(); 28 targetY[i] = m_vecAllVertex.at(i - vertexNumber)->pos.y(); 29 targetZ[i] = m_vecAllVertex.at(i - vertexNumber)->pos.z(); 30 } 31 else 32 { 33 //最高层的顶点坐标分量 34 targetX[i] = m_vecAllVertex.at((i - (vertexNumber + m_iAroundNumber) + (vertexNumber - m_iAroundNumber)))->pos.x(); 35 targetY[i] = m_vecAllVertex.at((i - (vertexNumber + m_iAroundNumber) + (vertexNumber - m_iAroundNumber)))->pos.y(); 36 targetZ[i] = m_vecAllVertex.at((i - (vertexNumber + m_iAroundNumber) + (vertexNumber - m_iAroundNumber)))->pos.z(); 37 } 38 } 39 40 Eigen::SimplicialCholesky<Eigen::SparseMatrix<float> > MatrixCholesky(MatrixLpls); 41 Eigen::VectorXf resX, resY, resZ; 42 resX = MatrixTranspose * targetX; 43 resY = MatrixTranspose * targetY; 44 resZ = MatrixTranspose * targetZ; 45 46 Eigen::VectorXf X, Y, Z; 47 X = MatrixCholesky.solve(resX); 48 //std::cout << X << std::endl; 49 Y = MatrixCholesky.solve(resY); 50 //std::cout << Y << std::endl; 51 Z = MatrixCholesky.solve(resZ); 52 //std::cout << Z << std::endl; 53 for (int i = 0; i < m_vecAllVertex.size(); i++) 54 { 55 pVERTEX pVertex = m_vecAllVertex.at(i); 56 float x, y, z; 57 x = X[i]; 58 y = Y[i]; 59 z = Z[i]; 60 61 pVertex->pos = osg::Vec3(X[i], Y[i], Z[i]); 62 } 63 }
参考文献：
Linear Rotation-invariant Coordinates for Meshes
Laplacian-Surface-Editing