图形处理（四）基于梯度场的网格编辑-Siggraph 2004

基于梯度场的网格编辑，对应的Paper为《Mesh Editing with Poisson-Based Gradient Field Manipulation》，是Siggraph 2004上的一篇paper，这篇paper与基于拉普拉斯的网格变形方法，统称为基于微分域的网格变形算法，这篇paper其实本质上最后的求解公式和基于拉普拉斯的网格变形方法一样，之所以能够siggraph，是因为它通过泊松梯度场的原理进行推导，算法的巧妙之处在于它以顶点(x,y,z)中的每一维作为一个标量场。

这篇paper涉及到的概念：散度、梯度场、标量场、向量场等看起来很难的东西，说实话，对于这篇paper因为网上找不到源代码，我把这篇paper看了好多遍，才把它的代码写出来。学这篇paper时是我第一次学习向量场的相关知识，向量场在三维算法中非常重要，同时当时给我的感觉也真不是一般的难，我看了好多关于标量场、矢量场的相关知识理论，才感觉慢慢理解。

一、相关理论

数学上的泊松方程：

其中f表示标量场，w表示梯度场。

三角网格曲面上的微分算子离散化(引用自《勾画式泊松网格编辑》)：

给定定义在网格曲面上的分段线性标量场f(v)=fi*φi(v),其中v为网格曲面上的任意一点；fi为标量场在网格曲面顶点vi处的函数值；φi(*)为分段线性基函数，它在顶点vi处取值为1 ，在其余顶点处取值为0。我们有标量场f 对应的梯度算子

其中▽φi(*)仅在顶点的邻接三角形上有非零值，且由于φi(*)分段线性，▽φi(*)在各个邻接三角形上为分段常值函数. 从几何角度，可以容易地给
出▽φi(*)在三角形T =(vi，vj ，vk)上的定义:

其中，R90代表绕三角形法向量nT 旋转90度，AT是三角形的面积。类似地，给定定义在三角网格曲面上的分段常值的矢量场w，我们定义在顶点vi处w的散度为：

根据梯度算子和散度算子的定义，最后可以推导出网格曲面上的标量场f在顶点vi处的拉普拉斯算子为：

这篇paper是由浙大的牛人周坤提出来的，算法最后跟拉普拉斯网格编辑的最后公式可以说是一样的，然而它的标量场给我很大的启示，这篇paper直接把(x,y,z)中的x，y，z分别当做一个标量场，然后对标量场求取梯度场，最后求取散度，然后通过泊松方程重建网格模型，实现网格变形。想要更深入的了解泊松重建，可以看看我的另外一篇博文《图像处理（十二）图像融合(1)Seamless cloning泊松克隆-Siggraph 2004》

二、算法实现

1、求取源网格曲面的梯度场，最后求取梯度场的散度。

[cpp] view plain copy

1. //计算各个顶点的梯度  
2. void CScaleDeformBrush::Get_Faces_Gradient()  
3. {  
4.     int fn=m_BaseMesh->faces.size();  
5.     m_BaseMesh->need_adjacentfaces();  
6.   
7.     #pragma omp parallel for  
8.     for (int i=0;i<fn;i++)  
9.     {  
10.         TriMesh::Face &f=m_BaseMesh->faces[i];  
11.         vec vij=m_BaseMesh->vertices[f[1]]-m_BaseMesh->vertices[f[0]];  
12.         vec vik=m_BaseMesh->vertices[f[2]]-m_BaseMesh->vertices[f[0]];  
13.         vec normalf=vij CROSS vik;  
14.         float areaf=0.5f*len(normalf);  
15.         normalize(normalf);  
16.         for (int k=0;k<3;k++)  
17.         {  
18.            m_Face_Gradient[i][k]=vec(0,0,0);  
19.            for (int j=0;j<3;j++)  
20.            {  
21.             vec ei=m_BaseMesh->vertices[f[(j+2)%3]]-m_BaseMesh->vertices[f[(j+1)%3]];  
22.             vec gradient=float(m_BaseMesh->vertices[f[j]][k]*0.5f/areaf)*(normalf CROSS ei);  
23.             m_Face_Gradient[i][k]=m_Face_Gradient[i][k]+gradient;  
24.             }  
25.   
26.         }  
27.     }  
28.   
29. }  
30. void CScaleDeformBrush::Compute_Divergence()  
31. {  
32.         //计算顶点的散度  
33.     m_BaseMesh->need_adjacentfaces();  
34.     int vn=m_BaseMesh->vertices.size();  
35.     #pragma omp parallel for  
36.     for (int i=0;i<vn;i++)  
37.     {     
38.         for (int j=0;j<3;j++)  
39.         {  
40.             m_vertices[i].VDivergence[j]=0.0f;  
41.         }  
42.         vector<int>&adjacentface=m_BaseMesh->adjacentfaces[i];  
43.         for (int j=0;j<adjacentface.size();j++)  
44.         {  
45.             TriMesh::Face &f=m_BaseMesh->faces[adjacentface[j]];  
46.             for (int k=0;k<3;k++)  
47.             {  
48.                 if (f[k]==i)  
49.                 {  
50.                     vec ei=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[(k+1)%3]];  
51.                     vec e1=m_BaseMesh->vertices[f[(k+1)%3]]-m_BaseMesh->vertices[f[k]];  
52.                     vec e2=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[k]];  
53.                     double cot_angle1=Cot_angle(e2,ei);  
54.                     double cot_angle2=Cot_angle(-1.0f*e1,ei);  
55.                     for (int xyz=0;xyz<3;xyz++)  
56.                     {  
57.                         m_vertices[i].VDivergence[xyz]+=0.5*(cot_angle1*(e1 DOT m_Face_Gradient[adjacentface[j]][xyz])+cot_angle2*(e2 DOT m_Face_Gradient[adjacentface[j]][xyz]));  
58.                     }  
59.                     break;  
60.                 }  
61.             }  
62.         }  
63.   
64.     }  
65.   
66. }  
67. //计算v1 v2 之间夹角的余切值  
68. double CScaleDeformBrush::Cot_angle(vec v1,vec v2)  
69. {  
70.     vec vivo=v1;  
71.     vec vjvo=v2;  
72.     double dotvector=vivo DOT vjvo;  
73.     dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);  
74.     return dotvector;  
75. }  

2、构建泊松方程的系数，矩阵A，也就是计算拉普拉斯矩阵

[cpp] view plain copy

1. //邻接顶点的余切权重计算  
2. void CScaleDeformBrush::CotangentWeights(TriMesh*TMesh,int vIndex,vector<double>&vweight,double &WeightSum,bool bNormalize)//计算一阶邻近点的各自cottan权重  
3. {     
4.     int NeighborNumber=TMesh->neighbors[vIndex].size();  
5.     vweight.resize(NeighborNumber);  
6.     WeightSum=0;  
7.     vector<int>&NeiV=TMesh->neighbors[vIndex];  
8.     for (int i=0;i<NeighborNumber;i++)  
9.     {  
10.         int j_nei=NeiV[i];  
11.         vector<int>tempnei;  
12.         Co_neighbor(TMesh,vIndex,j_nei,tempnei);  
13.         double cotsum=0.0;  
14.         for (int j=0;j<tempnei.size();j++)  
15.         {  
16.             vec vivo=TMesh->vertices[vIndex]-TMesh->vertices[tempnei[j]];  
17.             vec vjvo=TMesh->vertices[j_nei]-TMesh->vertices[tempnei[j]];  
18.             double dotvector=vivo DOT vjvo;  
19.             dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);  
20.             cotsum+=dotvector;  
21.         }  
22.         vweight[i]=cotsum/2.0;  
23.         WeightSum+=vweight[i];  
24.     }  
25.   
26.     if ( bNormalize )   
27.     {  
28.         for (int k=0;k<NeighborNumber;++k)  
29.         {  
30.             vweight[k]/=WeightSum;  
31.         }  
32.         WeightSum=1.0;  
33.     }  
34. }  
35.   
36. //获取两顶点的共同邻接顶点  
37. void CScaleDeformBrush::Co_neighbor(TriMesh *Tmesh,int u_id,int v_id,vector<int>&co_neiv)  
38. {  
39.     Tmesh->need_adjacentedges();  
40.     vector<int>&u_id_ae=Tmesh->adjancetedge[u_id];   
41.     int en=u_id_ae.size();  
42.     Tedge Co_Edge;  
43.     for (int i=0;i<en;i++)  
44.     {  
45.         Tedge &ae=Tmesh->m_edges[u_id_ae[i]];  
46.         int opsi=ae.opposite_vertex(u_id);  
47.         if (opsi==v_id)  
48.         {  
49.             Co_Edge=ae;  
50.             break;  
51.         }  
52.     }  
53.     for (int i=0;i<Co_Edge.m_adjacent_faces.size();i++)  
54.     {  
55.         TriMesh::Face af=Tmesh->faces[Co_Edge.m_adjacent_faces[i]];  
56.         for (int j=0;j<3;j++)  
57.         {  
58.             if((af[j]!=u_id)&&(af[j]!=v_id))  
59.             {  
60.                 co_neiv.push_back(af[j]);  
61.             }  
62.         }  
63.     }  
64. }  
65. //计算拉普拉斯矩阵  
66. void CScaleDeformBrush::Get_Laplace_Matrix()  
67. {  
68.     int vn=m_BaseMesh->vertices.size();  
69.     int count0=0;  
70.     vector<int>begin_N(vn);  
71.     for (int i=0;i<vn;i++)  
72.     {     
73.         begin_N[i]=count0;  
74.         count0+=m_BaseMesh->neighbors[i].size()+1;  
75.     }  
76.     typedef Eigen::Triplet<double> T;  
77.     std::vector<T> tripletList(count0);  
78.     for(int i=0;i<vn;i++)  
79.     {  
80.         VProperty & vi = m_vertices[i];  
81.         tripletList[begin_N[i]]=T(i,i,-vi.VSumWeight);  
82.         int nNbrs = vi.VNeighbors.size();  
83.         for (int k = 0;k<nNbrs;++k)   
84.         {  
85.             tripletList[begin_N[i]+k+1]=T(vi.VNeighbors[k],i,vi.VNeiWeight[k]);  
86.         }  
87.     }  
88.     m_Laplace_Matrix.resize(vn,vn);   
89.     m_Laplace_Matrix.setFromTriplets(tripletList.begin(), tripletList.end());  
90.   
91. }  

3、添加边界约束条件，并求解泊松方程，更新变形结果。

实时更新函数：

[cpp] view plain copy

1. void CScaleDeformBrush::Update_V_Position()  
2. {  
3.     Get_Faces_Gradient();//求梯度  
4.     int fn=m_BaseMesh->faces.size();  
5.     if(!m_ScaleFace.empty())  
6.     for (int i=0;i<fn;i++)  
7.     {  
8.         if(m_ScaleFace[i])  
9.         {  
10.             for (int j=0;j<3;j++)  
11.             {  
12.                 m_Face_Gradient[i][j]=1.1f*m_Face_Gradient[i][j];  
13.             }  
14.             m_BaseMesh->faces[i].beSelect=false;  
15.         }  
16.     }  
17.     Compute_Divergence();//求散度  
18.     if(!m_MatricesCholesky)//构建拉普拉斯矩阵  
19.     {  
20.         double a=m_Laplace_Matrix.coeff(0,0) +1;  
21.         m_Laplace_Matrix.coeffRef(0,0)=a;  
22.         m_MatricesCholesky=new Eigen::SimplicialCholesky<SparseMatrixType>(m_Laplace_Matrix);//矩阵分解  
23.     }  
24.     int vn=m_BaseMesh->vertices.size();  
25.     for (int i=0;i<3;i++)  
26.     {  
27.         Eigen::VectorXd rhs_xyz(vn);  
28.         for (int j=0;j<vn;j++)  
29.         {  
30.             rhs_xyz[j]=m_vertices[j].VDivergence[i];//方程组右边  
31.         }  
32.         rhs_xyz[0]=rhs_xyz[0]+1.0f*m_BaseMesh->vertices[0][i];  
33.         Eigen::VectorXd xyz=m_MatricesCholesky->solve(rhs_xyz);//求解方程  
34.         for (int j=0;j<vn;j++)  
35.         {  
36.             m_BaseMesh->vertices[j][i]=xyz[j];//更新结果  
37.         }  
38.     }  
39.     m_ScaleFace.clear();  
40.     m_ScaleFace.resize(fn,false);  
41.     m_BaseMesh->normals.clear();  
42.     m_BaseMesh->FaceNormal.clear();  
43.   
44. }  

接着我们来看一下用这个算法实现的简单局部编辑结果：

上面的实时局部缩放算法我是通过另外一篇paper《Differential-Based Geometry and Texture Editing with Brushes》的思想实现的，这篇paper基本上就是拷贝《Mesh Editing with Poisson-Based Gradient Field Manipulation》的思想，唯一的创新点在于它的实时交互设计方面，因为我是为了实现实时缩放刷，所以缩放的思想就是根据《Differential-Based Geometry and Texture Editing with Brushes》进行写代码的。

上面是用了上面的算法进行简单的实时编辑。

这篇paper后面还有后续的算法调整，比如梯度方向调整、还有实现网格融合、几何纹理Transfer。其中梯度方向调整是实现保特征变形的必备条件，因此如果你想要实现完整的算法，就要对梯度方向进行调整，这个可以参考我的以一篇博文《基于旋转不变量的网格变形》。在这里，我就不详细讲方向调整了，方向调整有专门的算法，paper很多。

须知：基于梯度域的变形方法和拉普拉斯网格变形算法一样，微分坐标不具有旋转不变的特点，在变形的时候，会发生曲面细节扭曲，需要对微分坐标，或者梯度方向进行调整，才能实现保特征变形，要实现旋转不变的变形，可以参考我的另外一篇博文《基于旋转不变量的网格变形》以此实现旋转不变的特点。

本文地址：http://blog.csdn.net/hjimce/article/details/46415291    作者：hjimce     联系qq：1393852684   更多资源请关注我的博客：http://blog.csdn.net/hjimce                  原创文章，转载请保留本行信息

参考文献：

1、《Mesh Editing with Poisson-Based Gradient Field Manipulation》

2、微分网格处理技术

3、勾画式泊松网格编辑