碼上歡樂
相關內容    簡體    繁體

梯度泊松方程網格形變

本文轉載自   查看原文   2018-02-23 20:04   900    圖形算法

　　網格頂點坐標的3個分量當做3個獨立的標量場，如此，三角面片便有3個獨立的梯度場，是為網格內在屬性。當頂點移動時，問題便為求解方程，

　　根據變分法得 其中Φ是為形變后的頂點坐標，W為形變后的梯度場。方程進一步用矩陣表示為 ，L為網格拉普拉斯矩陣，b為梯度場的散度.

　　1.定義梯度算子

　　設 f 為分段線性函數，網格三角面片{Xi, Xj, Xk}的頂點處有：f(Xi) = fi ; f(Xj) = fj ; f(Xk) = fk; u = (a, b)是三角面片的重心坐標，Bi, Bj, Bk是三角面片的一次拉格朗日

　　插值函數，通過線性插值，f 在三角面片上的梯度表示為。

　　函數Bi,Bj,Bk滿足Bi(u) + Bj(u) + Bk(u) = 1，所以Bi(u) + Bj(u) + Bk(u)的梯度等於0，所以 f 的梯度可表示為。

　　根據重心定理有，函數B的梯度為頂點到對邊的高的倒數，方向為對邊指向頂點。B的梯度表示為，A表示三角形面積，⊥表示旋轉90度。

　　2. 定義散度算子

　　設置向量值函數w：S——R3   。

　　S表示網格，W表示三角面片上的各個頂點向量，那么W的散度為。

　　3. 拉普拉斯算子

　　

　　

　　

　　

　構建拉普拉斯矩陣

  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         float totalWeight = 0.0;
  9         for (int j = 0; j < pVertex->vecNeighborEdge.size(); j++)
 10         {
 11             pEDGE pEdge = pVertex->vecNeighborEdge.at(j);
 12             pVERTEX pA = pEdge->pA;
 13             pVERTEX pB = pEdge->pB;
 14 
 15             pVERTEX pTarget;
 16             if (pA->id == pVertex->id)
 17                 pTarget = pB;
 18             else
 19                 pTarget = pA;
 20 
 21             std::vector<pTRIANGLE> vecTri = pEdge->vecNeighborTri;
 22             pVERTEX pC = NULL;
 23             float weight = 0.0;
 24             for (int k = 0; k < vecTri.size(); k++)
 25             {
 26                 pTRIANGLE pTri = vecTri.at(k);
 27                 pVERTEX p1 = pTri->pA;
 28                 pVERTEX p2 = pTri->pB;
 29                 pVERTEX p3 = pTri->pC;
 30                 if ((pA->id == p1->id && pB->id == p2->id) || (pA->id == p2->id && pB->id == p1->id))
 31                 {
 32                     pC = p3;
 33                 }
 34                 else if ((pA->id == p1->id && pB->id == p3->id) || (pA->id == p3->id && pB->id == p1->id))
 35                 {
 36                     pC = p2;
 37                 }
 38                 else if ((pA->id == p2->id && pB->id == p3->id) || (pA->id == p3->id && pB->id == p2->id))
 39                 {
 40                     pC = p1;
 41                 }
 42                 //開始求取權值
 43                 float cotAngle = 0.0;
 44                 GetCotAngle(pA->pos, pB->pos, pC->pos, cotAngle);
 45                 weight += cotAngle;
 46             }
 47             weight /= 2.0;
 48             pVertex->mapWeightForOther.insert(std::pair<int, float>(pTarget->id, weight) );
 49             totalWeight += weight;
 50         }
 51         pVertex->fTotalWeight = totalWeight;
 52     }
 53     //計算拉普拉斯坐標。方式二
 54     for (int i = 0; i < vertexNumber; i++)
 55     {
 56         pVERTEX pVertex = m_vecAllVertex.at(i);
 57         osg::Vec3 pos = pVertex->pos;
 58         pos = pos * pVertex->fTotalWeight;
 59         osg::Vec3 otherPos = osg::Vec3(0.0, 0.0, 0.0);
 60         for (int j = 0; j < pVertex->vecNeighborVertex.size(); j++)
 61         {
 62             pVERTEX pNeihbor = pVertex->vecNeighborVertex.at(j);
 63             std::map<int, float>::iterator itor = pVertex->mapWeightForOther.find(pNeihbor->id);
 64             float weight = itor->second;
 65             otherPos += pNeihbor->pos * weight;
 66         }
 67         pVertex->lplsPos = pos - otherPos;
 68     }
 69     
 70     int count = 0;
 71     std::vector<int> beginNumber(vertexNumber);
 72     for (int i = 0; i < vertexNumber; i++)
 73     {
 74         beginNumber[i] = count;
 75         count += m_vecAllVertex[i]->vecNeighborVertex.size() + 1;
 76     }
 77 
 78     std::vector<Eigen::Triplet<float> > vectorTriplet(count);
 79     for (int i = 0; i < vertexNumber; i++)
 80     {
 81         pVERTEX pVertex = m_vecAllVertex.at(i);
 82         //原始拉普拉斯矩陣
 83         vectorTriplet[beginNumber[i]] = Eigen::Triplet<float>(i, i, pVertex->fTotalWeight);
 84         //更改后
 85         //vectorTriplet[beginNumber[i]] = Eigen::Triplet<float>(i, i, -pVertex->fTotalWeight);
 86 
 87         int j = 0;
 88         std::map<int, float>::iterator itor;
 89         for(itor = pVertex->mapWeightForOther.begin(); itor != pVertex->mapWeightForOther.end(); itor++)
 90         {
 91             int neighborID = itor->first;
 92             float weight = itor->second;
 93             //原始拉普拉斯矩陣
 94             vectorTriplet[beginNumber[i] + j + 1] = Eigen::Triplet<float>(i, neighborID, -weight);
 95             //更改后
 96             //vectorTriplet[beginNumber[i] + j + 1] = Eigen::Triplet<float>(neighborID, i, weight);
 97             j++;
 98         }
 99     }
100 
101     //頭一圈和頂一圈
102     for (int i = 0; i < m_iAroundNumber * 2; i++)
103     {
104         float weight = 1.0;
105         pVERTEX pVertex;
106         if (i < m_iAroundNumber)
107         {
108             pVertex = m_vecAllVertex.at(i);
109         }
110         else
111         {
112             pVertex = m_vecAllVertex.at(i + m_iAroundNumber * 14);
113         }
114 
115         int row = i + vertexNumber;
116         vectorTriplet.push_back(Eigen::Triplet<float>(row, pVertex->id, weight));
117     }
118     
119     
120 
121     m_Matrix.resize(vertexNumber + m_iAroundNumber * 2, vertexNumber);
122     //m_Matrix.resize(vertexNumber, vertexNumber);
123     m_Matrix.setFromTriplets(vectorTriplet.begin(), vectorTriplet.end());
124     //std::cout << m_Matrix << std::endl;
125 }
View Code

　計算梯度

 1 void mainNode::CalculateGradient()
 2 {
 3     std::map<int, pTRIANGLE>::iterator itorOfTri;
 4     for (itorOfTri = m_mapAllTri.begin(); itorOfTri != m_mapAllTri.end(); itorOfTri++)
 5     {
 6         int triID = itorOfTri->first;
 7         pTRIANGLE pTri = itorOfTri->second;
 8         osg::Vec3 A = pTri->pA->pos;
 9         osg::Vec3 B = pTri->pB->pos;
10         osg::Vec3 C = pTri->pC->pos;
11 
12         
13         osg::Vec3 AB = B - A;
14         osg::Vec3 CA = A - C;
15         osg::Vec3 BC = C - B;
16 
17         osg::Vec3 normal = ((B - A) ^ (C - A));
18         //三角形面積
19         float area = normal.length() * 0.5;
20         normal.normalize();
21 
22         float times = 2.0;
23         osg::Vec3 vecGradientX = ((normal ^ BC) / (times * area) * A.x()) + ((normal ^ AB) / (times * area) * C.x()) + ((normal ^ CA) / (times * area) * B.x());
24         osg::Vec3 vecGradientY = ((normal ^ BC) / (times * area) * A.y()) + ((normal ^ AB) / (times * area) * C.y()) + ((normal ^ CA) / (times * area) * B.y());
25         osg::Vec3 vecGradientZ = ((normal ^ BC) / (times * area) * A.z()) + ((normal ^ AB) / (times * area) * C.z()) + ((normal ^ CA) / (times * area) * B.z());
26         pTri->vecGradient[0] = vecGradientX;
27         pTri->vecGradient[1] = vecGradientY;
28         pTri->vecGradient[2] = vecGradientZ;
29     }
30 }
View Code

　　計算散度

 1 void mainNode::CalculateDivergence()
 2 {
 3     
 4     for (int i = 0; i < m_vecAllVertex.size(); i++)
 5     {
 6         pVERTEX pVertex = m_vecAllVertex.at(i);
 7         std::vector<pTRIANGLE> vecTri= pVertex->vecNeighborTri;
 8 
 9         float divergenceX, divergenceY, divergenceZ;
10         divergenceX = divergenceY = divergenceZ = 0.0;
11 
12         for (int j = 0; j < vecTri.size(); j++)
13         {
14             pTRIANGLE pTri = vecTri.at(j);
15             pVERTEX pA = pTri->pA;
16             pVERTEX pB = pTri->pB;
17             pVERTEX pC = pTri->pC;
18             //A為主頂點
19             osg::Vec3 A, B, C;
20             if (pA->id == pVertex->id)
21             {
22                 A = pA->pos;
23                 B = pB->pos;
24                 C = pC->pos;
25             }
26             else if (pB->id == pVertex->id)
27             {
28                 A = pB->pos;
29                 B = pC->pos;
30                 C = pA->pos;
31             }
32             else if (pC->id == pVertex->id)
33             {
34                 A = pC->pos;
35                 B = pA->pos;
36                 C = pB->pos;
37             }
38             //
39             osg::Vec3 AB, AC;
40             float cotAngleABC, cotAngleACB;
41             AB = B - A;
42             AC = C - A;
43             GetCotAngle(A, C, B, cotAngleABC);
44             GetCotAngle(A, B, C, cotAngleACB);
45 
46             for (int k = 0; k < 3; k++)
47             {
48                 float divergenceTemp = 0 - 0.5 * (cotAngleABC * (AC * pTri->vecGradient.at(k)) + cotAngleACB * (AB * pTri->vecGradient.at(k)));
49                 if (k == 0)
50                 {
51                     divergenceX += divergenceTemp;
52                 }
53                 else if (k == 1)
54                 {
55                     divergenceY += divergenceTemp;
56                 }
57                 else if(k == 2)
58                 {
59                     divergenceZ += divergenceTemp;
60                 }
61             }
62         }
63 
64         pVertex->vecDivergence[0] = divergenceX;
65         pVertex->vecDivergence[1] = divergenceY;
66         pVertex->vecDivergence[2] = divergenceZ;
67     }
68 
69 }
View Code

 

　　參考論文：

　　Mesh Editing with Poisson-Based Gradient Field Manipulation

　　三角網格曲面上的Laplace算子及其應用

　　http://blog.sina.com.cn/s/blog_73428e9a0101h4l9.html

　

 

　　

　　

　　


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 泊松方程解法 從泊松方程的解法，聊到泊松圖像融合 泊松分布、泊松過程、泊松點過程 泊松分布與泊松回歸模型 泊松分布 泊松融合 泊松分布 圖像的泊松(Poisson)編輯、泊松融合 Poisson Distribution——泊松分布 泊松分布到底是什么
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM