基于梯度场的网格编辑，对应的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 plaincopy
  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，也就是计算拉普拉斯矩阵

[cpp] view plaincopy
  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、添加边界约束条件，并求解泊松方程，更新变形结果。

实时更新函数：

[cpp] view plaincopy
  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很多。

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

本文地址：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、勾画式泊松网格编辑

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

  1. 图形处理（三）简单拉普拉斯网格变形-Siggraph 2004

    三角网格变形一直是CAGD相关领域的重点,刚上研究生的时候,感觉有点神奇.而且一上来导师就给我发了一篇基于格林坐标的自由变形的相关paper,让我看,外文文献,看了n多天,第一次看外文文献,啥也没看懂 ...

  2. 图形处理（五）基于旋转不变量的网格变形-Siggraph 2007

    一.相关理论 本篇博文主要讲解2007年Siggraph上的一篇经典paper:<Linear Rotation-invariant Coordinates for Meshes>,这篇p ...

  3. 图形处理（十二）拉普拉斯网格优化、最小二乘网格模型光顺

    看这篇博文前,请先参考我的另外一篇博文<图形处理(三)简单拉普拉斯网格变形-Siggraph 2004>学习拉普拉斯坐标的相关理论知识.这里要分享的paper,是通过拉普拉斯的方法实现三角 ...

  4. 机器学习笔记之玻尔兹曼机(三)梯度求解(基于平均场理论的变分推断)

    机器学习笔记之玻尔兹曼机--基于平均场推断梯度求解 引言 回顾:玻尔兹曼机模型参数梯度求解困难与MCMC方法的处理方式 变分推断方法处理玻尔兹曼机对数似然梯度 引言 上一节介绍了使用马尔可夫链蒙特卡洛 ...

  5. 机器学习笔记之玻尔兹曼机(三)基于平均场理论变分推断的梯度求解(续)

    机器学习笔记之玻尔兹曼机--基于平均场推断梯度求解[续] 引言 Λ 3 \Lambda_3 Λ3​梯度求解 求解最优参数 ϕ ^ j \hat {\phi}_j ϕ^​j​ 引言 基于玻尔兹曼机(三) ...

  6. 吴恩达机器学习（十四）推荐系统（基于梯度下降的协同过滤算法）

    目录 0. 前言 1. 基于内容的推荐算法(Content-based recommendations) 2. 计算电影特征 3. 基于梯度下降的协同过滤算法(Collaborative filter ...

  7. 图形处理（七）基于热传播的测地距离计算-Siggraph 2013

    这里要跟大家分享的是2013年Siggraph上面的一篇paper,名为<Geodesics in Heat:A New Approach to Computing Distance Based ...

  8. 【数据库实验】实验四 基于嵌入SQL的综合应用编程（基于QSqlTableModel实现）

    [数据库实验]实验四 基于嵌入SQL的综合应用编程 一.实验目的 二.实验要求 三.实验内容.实验结果与主要程序代码 数据准备(建表并插入数据) 前言:黎哥的写法 参考 建表 SQL语言插入数据 S表 ...

  9. 云原生爱好者周刊：Linkerd 即将赢得这场服务网格战争的胜利？

    云原生一周动态要闻: 供应链安全项目 in-toto 升级成为 CNCF 孵化项目 云原生可观测性微调查结果出炉 Linkerd 发布 Kubernetes 自动多集群故障转移新特性 Solo.io ...

最新文章

  1. mac 拷贝文件时报错 8060 解决方案
  2. php windows 编译,Windows编译PHP7.2拓展
  3. Testin云测试：QQ(4.2.0)安卓版客户端可用性优秀
  4. python可以处理什么文件夹_Python处理文件和文件夹的10条命令
  5. 混凝土墙开洞_易县混凝土剪力墙切割常见问题
  6. 数据采集时总提示未登录_个税申报系统新功能！申报数据丢了也能找回！|税务局|个税|办税服务厅|纳税...
  7. SQL Server 2008 R2 SSRS 安装配置后无法使用问题的解决方法
  8. 使用ssl_exporter监控K8S集群证书
  9. 为枪击事件默哀，程序员们确实要重视代码规范
  10. 【IIOT】欧姆龙PLC数采之NX/NJ系列
  11. 使用Photoshop制作名片
  12. RADIUS协议指南
  13. thinkphp 后台管理框架swiftadmin的使用
  14. ps抠图怎么放大图片_PS抠图时选区图片放大后,怎么移动图片抠图选区？
  15. C51最小单片机系统
  16. 基于IPv6的5G专网终端身份认证技术与应用
  17. 曲率、曲率（对弧长）的导数以及曲率导数（对弧长）的导数的计算
  18. Caused by: java.util.concurrent.ExecutionException: java.util.concurrent.ExecutionException: com.and
  19. [MATLAB] 心形图
  20. DFS和BFS求字符串的所有非空子集———Java

热门文章

  1. 2014年应该学习的十种编程语言
  2. 食品安全溯源区块链解决方案探索-转载
  3. @excel 注解_7 行代码实现 Excel 文件导出
  4. Shell-实际业务操作02
  5. Oracle优化12-10053事件
  6. CoordinatorLayout+AppBarLayout实现上滑隐藏ToolBar-Android M新控件
  7. linux接收网络数据并存存储,linux网络数据包数据结构 Socket Buffer
  8. java shape用法_Java PShape.scale方法代码示例
  9. 一个对象的属性_【前端冷知识】如何判断一个对象的某个属性是可写的？
  10. 汇编编程计算机流程图,汇编程序怎么做流程图？