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 的旋轉矩陣 Tij

　　關于矩陣 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 } View Code

　　構建旋轉矩陣方程組

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 } View Code

　　求取新的局部標架并獲得新的拉普拉斯坐標

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 } View Code

　　構建拉普拉斯矩陣

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 } View Code

　　根據新的拉普拉斯坐標求取新的頂點坐標

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 } View Code

?

　　參考文獻：

　　Linear Rotation-invariant Coordinates for Meshes

　　Laplacian-Surface-Editing

　　http://blog.csdn.net/hjimce/article/details/46415435

轉載于:https://www.cnblogs.com/TooManyWayToBe/p/8465470.html

總結

以上是生活随笔為你收集整理的协同旋转不变网格形变的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。