补洞与fairing算法

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

学习文章 关注公众号回复 Filling Holes in Meshes可获取电子书

补洞算法

作用

• 对空洞区域（洞，缝细，多个孤岛）进行网格填补
• 优化填补的网格三角形
• 平滑网格，达到与周边已有三角形1，2，3阶平滑连续

整体步骤

1. 识别洞
2. 洞三角化
3. 重新网格化，优化三角形
4. 洞光顺

识别洞

基本原理是遍历每个点，如果发现是边界点，则搜索与其连通的所有边界点，构成一个洞。

hole bfs(mesh, v) {
	hole = [];
	que.push(v);
	vis[v]=true;
	while(!que.empty()){
		f = que.front();
		que.pop();
		hole.add(v);
		for each vertex v1 in 1ring(f) {
			if(!vis[v1] && is_boundary(v1)){
				que.push(v1);
				vis[v1]=true;
			}
		}
	}
	return hole;
}

holes detect_holes(mesh) {
	holes = [];
	for each v inmesh {
		if(!vis[v] && is_boundary(v)) holes.add(get_hole(v));	
	}
	return holes;
}

动态规划三角化

一个洞可以认为是一堆逆时针排列的点组成的立体多边形p(v0, v1, vi,…vj,…vn-1)。

上图一共有n个点。
我们设dp[i][j]为以第i点为起点逆时针到第j点为终点的多边形三角化后的权值结果，
那么dp[0][n-1]就是是整体权值结果，我们要去最优化dp[0][n-1]。

对于多边形Q(vi,vi+1, …vj)来说，我们知道边（i,j）最终肯定要成为一个三角的边。

那么我们就枚举i<m<j，剩下的点就被分成了R(vi, vi+1, …vm), L(vm, vm+1, …vj).

d p [ i ] [ j ] = m i n i < m < j { d p [ i ] [ m ] + d p [ m ] [ j ] + ϕ ( v i , v m , v j ) } , ϕ 为权值计算函数，有多种策略可以选择 dp[i][j] = \underset{i<m<j}{min} \{dp[i][m]+dp[m][j] +\phi(vi, vm,vj) \}, \phi为权值计算函数，有多种策略可以选择 dp[i][j]=i<m<jmin​{dp[i][m]+dp[m][j]+ϕ(vi,vm,vj)},ϕ为权值计算函数，有多种策略可以选择

权值策略

• 面积
  $\varphi (vi,vm,vj)为三点组成三角形的面积$
  那么dp[0][n-1] 就是取面积最小

• 二面角+面积
  $\varphi (vi,vm,vj)为(a,b)=（最大二面角，面积）二值对$
  那么dp[0][n-1] 就是（最大二面角最小，面积最小）
  二值对比较规则
  (a, b) <(c, d) 当a<c || (a==c && b<d)
  二值对运算规则
  (a, b) + (c, d) = (max(a, c), b+d)
  二面角计算规则
  
  $上图中\varphi (vi,vm,vj)的二面角值为2与1的二面角和3与1的二面角最大值$

动态规划实现

void dp_triangulation(polygon p) {
	select[][]// 存储select[i][j]要选择哪个点来连接
	dp[][] // 存储最佳权值
	for i=0 to n {dp[i][i+1]=0}
	
	// 枚举长度从小到大
	for len=2 to n {
		for i=0 to n {
			j=i+len
			for m=i+1 to j {
				if(dp[i][m] + dp[m][j] + phi(i,m,j)<dp[i][j]) {
					dp[i][j] = dp[i][m] + dp[m][j] + phi(i,m,j)
					select[i][j] =m
				}
			}
		}
	}
	
	// 添加边
	q.push({0, n-1}) // 添加区域
	while(!q.empty()) {
		f = q.front();
		q.pop();
		m = select[f.first][f.second]
		if(m>=0) {
			p.addedge(f.first, m)
			p.addedge(m, f.second)
			// 分列区间
			q.push({f.first, m})
			q.push({m, f.second})
		}
	}
}

重新网格化，优化三角形

relex 操作

对于1条边的两个三角形，如果（不满足空圆性）某个三角形的外接球内有相邻三角形的点，则需要flip边。

1.对于所有洞边界的点，计算平均边长avl[v]
2.遍历每个三角形, 点为i,j,k, 中点为vc, avl[vc] = (avl[vi]+avl[vj]+avl[vk])/3
3. 如果对于m对于任意i,j,k 满足a|vc-vm|>alv[vc] && a|vc-vm|>alv[vm]，则分裂三角，并relax新加的三条边。
4. 对所有内部点relax
5. 循环4，2，直到网格不发生变化

fairing

到上一步为止我们在洞区域可以得到一个比较稠密且均匀的三角形，但是形状上还不光滑，和旁边的几何形状也格格不入。

fairing的目的就是把洞区域平滑，还可以与旁边的几何形状有很好的过度。

雨伞算子

令w(vi,vj) 为 边（vi, vj）的权值，w(vi)为vi1邻域所有权值之和。

可以定义一个

带权重的雨伞算子（形状如雨伞，红点为伞顶v）

$U(v)=-v+\frac{1}{w(v)}\sum _{i}w(v,v_i)v_i,w(v)=\sum _{i}w(v,v_i)$

平滑的常规做法就是让所有U(v)=0, 对于洞边界的点尽量保持原来的值，组成Av=b的形式，来优化v。

权重计算策略的不同会带来不同的效果，

如果 所有w(vi,vj)=1， 则称为uniform 雨伞算子，会得到一个平均曲面。
如果 w(vi,vj)= |vi-vj|， 则称为 缩放雨伞算子，会得到一个缩放效果。
如果


则称为harmonic 雨伞算子。

除了一阶算子，还有二阶算子，可以让曲面过度得更平滑

二阶算子就是对1阶算子做一次加权

$U^2(v)=-U(v)+\frac{1}{w(v)}\sum _{i}w(v,v_i)U(v_i),w(v)=\sum _{i}w(v,v_i)$

不同的平滑效果对应如下

上图分别为一阶，二阶，三阶平滑效果

方程推导

一阶平滑

我们要求对于洞边界处以及至少1邻域的点（原来网格上有的点），保持不动。

蓝色部分不动，黄色部分为fairing区域(U(v)=0)，为了方便计算，我将整个网格都参与平滑

假设整个洞相关的点有n个（包含蓝色部分），其中蓝色部分为前m个点
则与蓝色部分相关的有2条方程（一条为固定方程在A的最后m行，一条为平滑方程）
黄色部分有1条方程（平滑方程）
总方程数为n+m条。

Av=b 形式如下

$\left[\begin{array}{cccc}-1& \frac{w(v1,v2)}{\sum _{i}w(v_1,v_i)}& ...& \frac{w(v1,vn)}{\sum _{i}w(v_1,v_i)}\\ \frac{w(v2,v1)}{\sum _{i}w(v_2,v_i)}& -1& ...& \frac{w(v2,vn)}{\sum _{i}w(v_2,v_i)}\\ ...& ...& ...& ...\\ \frac{w(vn,v1)}{\sum _{i}w(v_n,v_i)}& ...& \frac{w(vn,vn-1)}{\sum _{i}w(v_n,v_i)}& -1\\ 1& 0& ...& 0\\ 0& 1& ...& 0\\ 0& 0& 1...& 0\\ ...& ...& ...& 0\end{array}\right]{\left[\begin{array}{c}{v}_1^{\prime }\\ v_2^{\prime }\\ ...\\ v_n^{\prime }\end{array}\right]}^T=\left[\begin{array}{c}0\\ 0\\ ...\\ 0\\ v_1\\ v_2\\ ...\\ v_m\end{array}\right]$

上式中A的前n行为平滑优化方程，后面m行的前m列为单位矩阵，表示希望前m个点移动得越少越好。

优化上述方程就可以得到平滑效果。

二阶平滑

原理与一阶平滑类似，A的后m行和之前一样，只是把U(v)=0替换成U^2(v)=0。

$将U^2(v)展开得到：$

$U^2(v_i)=v_i-\frac{1}{w(v_i)}\sum _{j}w(v_i,v_j)v_j+\frac{1}{w(v_i)}\sum _{j}w(v_i,v_j)\left[-v_j+\frac{1}{w(v_j)}\sum _{k}w(v_j,v_k)v_k\right]$

$=v_i+\frac{1}{w(v_i)}\sum _{j}w(v_i,v_j)\left[-2v_j+\frac{1}{w(v_j)}\sum _{k}w(v_j,v_k)v_k\right]$

根据公式写出A矩阵前n行构造的伪代码

matrix gen_matrix() {
	matrix A
	for i=0 to n {
		A[i][i]+=1
		for each vi's onering vj {
			w1 = w(vi, vj)/w(vi)
			A[i][j]-=2w1
			for each vj's onering vk {
				w2 = w(vj, vk)/w(vj)
				A[i][k] += w1*w2
			}
		}
	}
}

注意点：由于上述都是优化性质的，所以不能保证原来的固定点不动。可以根据实际需要将原来的点拉回来。

---

本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。创作不易，帮忙点击公众号的链接。