计算机视觉
图像处理

图形处理(四)基于梯度场的网格编辑-Siggraph

基于梯度场网格编辑,对应的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、求取源网格曲面的梯度场,最后求取梯度场的散度。

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

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

  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.     if ( bNormalize )
  26.     {
  27.         for (int k=0;k<NeighborNumber;++k)
  28.         {
  29.             vweight[k]/=WeightSum;
  30.         }
  31.         WeightSum=1.0;
  32.     }
  33. }
  34. //获取两顶点的共同邻接顶点
  35. void CScaleDeformBrush::Co_neighbor(TriMesh *Tmesh,int u_id,int v_id,vector<int>&co_neiv)
  36. {
  37.     Tmesh->need_adjacentedges();
  38.     vector<int>&u_id_ae=Tmesh->adjancetedge[u_id];
  39.     int en=u_id_ae.size();
  40.     Tedge Co_Edge;
  41.     for (int i=0;i<en;i++)
  42.     {
  43.         Tedge &ae=Tmesh->m_edges[u_id_ae[i]];
  44.         int opsi=ae.opposite_vertex(u_id);
  45.         if (opsi==v_id)
  46.         {
  47.             Co_Edge=ae;
  48.             break;
  49.         }
  50.     }
  51.     for (int i=0;i<Co_Edge.m_adjacent_faces.size();i++)
  52.     {
  53.         TriMesh::Face af=Tmesh->faces[Co_Edge.m_adjacent_faces[i]];
  54.         for (int j=0;j<3;j++)
  55.         {
  56.             if((af[j]!=u_id)&&(af[j]!=v_id))
  57.             {
  58.                 co_neiv.push_back(af[j]);
  59.             }
  60.         }
  61.     }
  62. }
  63. //计算拉普拉斯矩阵
  64. void CScaleDeformBrush::Get_Laplace_Matrix()
  65. {
  66.     int vn=m_BaseMesh->vertices.size();
  67.     int count0=0;
  68.     vector<int>begin_N(vn);
  69.     for (int i=0;i<vn;i++)
  70.     {
  71.         begin_N[i]=count0;
  72.         count0+=m_BaseMesh->neighbors[i].size()+1;
  73.     }
  74.     typedef Eigen::Triplet<double> T;
  75.     std::vector<T> tripletList(count0);
  76.     for(int i=0;i<vn;i++)
  77.     {
  78.         VProperty & vi = m_vertices[i];
  79.         tripletList[begin_N[i]]=T(i,i,-vi.VSumWeight);
  80.         int nNbrs = vi.VNeighbors.size();
  81.         for (int k = 0;k<nNbrs;++k)
  82.         {
  83.             tripletList[begin_N[i]+k+1]=T(vi.VNeighbors[k],i,vi.VNeiWeight[k]);
  84.         }
  85.     }
  86.     m_Laplace_Matrix.resize(vn,vn);
  87.     m_Laplace_Matrix.setFromTriplets(tripletList.begin(), tripletList.end());
  88. }

3、添加边界约束条件,并求解泊松方程,更新变形结果。实时更新函数:

  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. }

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

上面的实时局部缩放算法我是通过另外一篇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很多。

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

参考文献:

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

2、微分网格处理技术

3、勾画式泊松网格编辑

转载注明来源:CV视觉网 » 图形处理(四)基于梯度场的网格编辑-Siggraph

分享到:更多 ()
扫描二维码,给作者 打赏
pay_weixinpay_weixin

请选择你看完该文章的感受:

0不错 0超赞 0无聊 0扯淡 0不解 0路过

评论 2

评论前必须登录!