欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

本文实现的算法来自下面的论文
参考文献Bilateral Normal Filtering for Mesh Denoising ，https://ieeexplore.ieee.org/document/5674028

算法原理及过程

算法原理

先利用双边滤波对面法向进行调整，然后根据法向去调整顶点，使得面的法向尽量接近调整后的法向。

法向调整

$n^{k+1}(f_i)=\frac{1}{K_p}\sum _{j\in \Omega (i)}{A}_j{W}_s(||c_i-c_j||)W_r(||n^k(f_i)-n^k(f_j)||)\cdot n^k(f_j)$

上述式子中

k代表迭代次数

$n^k(f_i)表示经过k次迭代后第i个面的法向$

$c_i代表第i个面的中心$

$\Omega (i)表示第i个面的所有邻面$

$A_j表示第j个面的面积$

$$\begin{gathered} W_s(x)是一个高斯函数=e^{-\frac{x^2}{2{\sigma }_s^2}}, \\ 其中x=c_i-c_j的长度，{\sigma }_s是可调参数， \\ 本次实现中取值为所有邻面中心差值的平均值 \end{gathered}$$

$$\begin{gathered} W_r(x)是一个高斯函数=e^{-\frac{x^2}{2{\sigma }_r^2}}, \\ 其中x=n^k(f_i)-n^k(f_j)的长度，{\sigma }_r是可调参数， \\ 本次实现中取值为1 \end{gathered}$$

$K_p=\sum _{j\in \Omega (i)}{A}_j{W}_s(||c_i-c_j||)W_r(||n^k(f_i)-n^k(f_j)||)$

法向调整可以运用多次，具体参数可调k。
每次计算完法向后，要对法向进行归一化。

点顶恢复

想要根据法向计算出所有顶点是非常困难的。这里采用单步高斯迭代的方式，每次只调整一个点顶，使得法向尽可能接近计算出来的法向。

根据定义我们有如下公式

$\{\begin{array}{l}{n}_f\cdot (x_j-x_i)=0\\ n_f\cdot (x_k-x_j)=0\\ n_f\cdot (x_i-x_k)=0\end{array}$

$x_i,x_j,x_k是某个面的三点，n_f是调整后的法向，上述公式表示每条边都要垂直于法线$

$我们的目标是使上述公式越接近于0越好$

可以得到一个能量方程

$E=\sum _{f_k}\sum _{i,j\in f_k}(n_k\cdot (x_j-x_i){)}^2$

现在要做的就是求一组顶点，使得E最小。

这里每次只移动一个点，使得E最小。

对上述方程进行求导并使其等于0，可以推导出某一点的计算公式

$x_i^{new}=x_i+\frac{1}{N_i}\sum _{f_j\in \Omega (i)}{n}_j\cdot (n_j\cdot (c_j-x_i))$

$N_i表示顶点的度$

$\Omega (i)是第i个点所有邻面$

$x_i是原来要变化的点$

$c_j是其他顶点的坐标$

$n_j是第j个面的目标法向$

上述过程每次对所有顶点进行恢复，可以进行多次参数可调，本次实现采用8次。

算法实现

代码github https://github.com/LightningBilly/ACMAlgorithms/tree/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw3

#include <iostream>
#include <fstream>
#include <sstream>
#include "../../PolyMesh/IOManager.h"
#include <string>
#include <unistd.h>using namespace std;
using namespace acamcad;
using namespace polymesh;void hw3()
{puts("hw3");char buffer[500];getcwd(buffer, 500);printf("The current directory is: %s/../\n", buffer);string mesh_path = buffer;mesh_path += "/../src/hw3/bunny_random.obj";double SigmaCenter, SigmaNormal=1;double NormalIterNum=8, VertexIterNum=8;PolyMesh* mesh = new PolyMesh();loadMesh(mesh_path, mesh);mesh->updateMeshNormal();std::vector<MVector3> NewNormal(mesh->numPolygons());std::vector<double> FaceArea(mesh->numPolygons());std::vector<MPoint3> FaceCenter(mesh->numPolygons());for (MPolyFace* fh : mesh->polyfaces()){int f_id = (*fh).index();NewNormal[f_id] = (*fh).normal();std::vector<MVert*> P;for (FaceVertexIter vv_it = mesh->fv_iter(fh); vv_it.isValid(); ++vv_it){P.push_back(*vv_it);}auto e12 = P[1]->position() - P[0]->position();auto e13 = P[2]->position() - P[0]->position();double area = cross(e12, e13).norm() * 0.5;FaceArea[f_id] = area;FaceCenter[f_id] = mesh->calculatFaceCenter(fh);}/** 计算点标准差参数* 把所有邻面的中点作差* 最后除以面的3倍*/SigmaCenter = 0;for (MPolyFace* fh : mesh->polyfaces()){int f_id = (*fh).index();for (FaceFaceIter nei_fh = mesh->ff_iter(fh); nei_fh.isValid(); ++nei_fh){int ff_id = (*nei_fh)->index();SigmaCenter += (FaceCenter[f_id] - FaceCenter[ff_id]).norm();}}SigmaCenter /= mesh->numPolygons() * 3;/** 利用法向迭代公式计算法向*/for (int i = 0; i < NormalIterNum; i++){for (MPolyFace* fh : mesh->polyfaces()){double Kp = 0;MVector3 NewN(0, 0, 0);int fh_id = (*fh).index();for (FaceFaceIter nei_fh = mesh->ff_iter(fh); nei_fh.isValid(); ++nei_fh){int nei_fh_id = (*nei_fh)->index();double delta_center = (FaceCenter[fh_id] - FaceCenter[nei_fh_id]).norm();double delta_normal = (NewNormal[fh_id] - NewNormal[nei_fh_id]).norm();double Aj = FaceArea[nei_fh_id];double Ws = exp(-delta_center * delta_center / (2 * SigmaCenter * SigmaCenter));double Wr = exp(-delta_normal * delta_normal / (2 * SigmaNormal * SigmaNormal));NewN += Aj * Ws * Wr * NewNormal[nei_fh_id];Kp += Aj * Ws * Wr;}NewNormal[fh_id] = NewN / Kp;NewNormal[fh_id] /= NewNormal[fh_id].norm();}}// 调整点位置for (int i = 0; i < VertexIterNum; i++){for (MVert* vh : mesh->vertices()){MPoint3 x_i = (*vh).position();MPoint3 delta_xi(0, 0, 0);int Nei_count = 0;for (VertexFaceIter fh = mesh->vf_iter(vh); fh.isValid(); ++fh){Nei_count++;MPoint3 cj = mesh->calculatFaceCenter(*fh);MVector3 nj = NewNormal[(*fh)->index()];delta_xi = delta_xi + nj * (nj.data()[0] * (cj - x_i).data()[0] + nj.data()[1] * (cj - x_i).data()[1] + nj.data()[2] * (cj - x_i).data()[2]);}x_i = x_i + delta_xi / Nei_count;(*vh).setPosition(x_i);}}// 输出网格数据writeMesh("result.obj", mesh);
}

算法效果

原型网格 https://github.com/LightningBilly/ACMAlgorithms/blob/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw3/bunny_random.obj

平滑后的网格
https://github.com/LightningBilly/ACMAlgorithms/blob/master/%E5%9B%BE%E5%BD%A2%E5%AD%A6%E7%AE%97%E6%B3%95/%E4%B8%89%E8%A7%92%E7%BD%91%E6%A0%BC%E7%AE%97%E6%B3%95/MeshWork/src/hw3/result.obj

原模型

平滑后模型

从上图中可以看到有很明显平滑效果。

---

本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。