使用Quadric Error Metrics进行mesh简化

背景

Hoppe等人在1993年的时候首先提出了一系列的mesh optimization技术[6],其中包括edge collapse, edge split, edge swap等操作。基于edge的操作相比较于之前的模型简化(或者说down-sampling)技术,如triangel collapse,最大的好处是能保证模型原始的拓扑结构,因此可以进一步用于其他的geometry processing操作。
关于模型化简时选择最合适的边进行collapse,我们需要选择那些坍缩后的边对原有的模型影响最小的边。Hoppe在论文中,使用的是能量方程来判断。但实际情况中,计算的成本过高。因此,在1997年的时候,Garland和Heckbert等人提出了一种基于Quaderic Error的度量标准[5]来判断最合适的顶点。

理论推导

点距平面的误差距离

我们用一个顶点到它的1-ring相邻面的距离之和来表示它的误差。

点v到平面上一点p的距离如图所示

我们采用齐次坐标系下的向量形式表示顶点v
$$
v=\begin{bmatrix}
v_{x}\\
v_{y}\\
v_{z}\\
1
\end{bmatrix}
$$
同时,平面表达式$ax+by+cy+d=0$写成矩阵形式就是
$$
\begin{array}{l}
\begin{bmatrix}
a & b & c & d
\end{bmatrix}\begin{bmatrix}
p_{x}\\
p_{y}\\
p_{z}\\
1
\end{bmatrix} =0\\
\Longrightarrow \begin{bmatrix}
a & b & c & 0
\end{bmatrix}\begin{bmatrix}
p_{x}\\
p_{y}\\
p_{z}\\
1
\end{bmatrix} =-d
\end{array}
$$
前面提到的距离表达式可以写成
$dist( v) \ =\ dot( N,\ v-p) =dot( N,v) -dot( N,p)$
平面法线由于是方向向量,故在齐次坐标系下w分量为0:$\begin{bmatrix}a & b & c & 0\end{bmatrix}^{T}$
由前面的几个式子,从而可以进一步化简$dist(v)$
$$
\begin{array}{l}
=\begin{bmatrix}
a & b & c & 0
\end{bmatrix}\begin{bmatrix}
v_{x}\\
v_{y}\\
v_{z}\\
1
\end{bmatrix} -\begin{bmatrix}
a & b & c & 0
\end{bmatrix}\begin{bmatrix}
p_{x}\\
p_{y}\\
p_{z}\\
1
\end{bmatrix}\\
=\begin{bmatrix}
a & b & c & 0
\end{bmatrix}\begin{bmatrix}
v_{x}\\
v_{y}\\
v_{z}\\
1
\end{bmatrix} +d\\
=\begin{bmatrix}
a & b & c & d
\end{bmatrix}\begin{bmatrix}
v_{x}\\
v_{y}\\
v_{z}\\
1
\end{bmatrix}\\
\Longleftrightarrow \begin{bmatrix}
v_{x} & v_{y} & v_{z} & 1
\end{bmatrix}\begin{bmatrix}
a\\
b\\
c\\
d
\end{bmatrix}
\end{array}
$$
$dist(v)$可正可负,我们只关注绝对值,所以我们对它求平方
$$
dist( v)^{2} =\left(\begin{bmatrix}
v_{x} & v_{y} & v_{z} & 1
\end{bmatrix}\begin{bmatrix}
a\\
b\\
c\\
d
\end{bmatrix}\right) .\left(\begin{bmatrix}
a & b & c & d
\end{bmatrix}\begin{bmatrix}
v_{x}\\
v_{y}\\
v_{z}\\
1
\end{bmatrix}\right)
$$
利用矩阵的结合律,$dist( v)^{2}$ 可以进一步表示为
$$
=v^{T}\begin{bmatrix}
a\\
b\\
c\\
d
\end{bmatrix}\begin{bmatrix}
a & b & c & d
\end{bmatrix} v
$$

$$
Q_p=\begin{bmatrix}
a\\
b\\
c\\
d
\end{bmatrix}\begin{bmatrix}
a & b & c & d
\end{bmatrix}
$$
这是一个只和平面相关的参数,可以表示任意一点到该平面的距离系数。我们这里将其称为点v到平面P的误差矩阵
$$
Q_{p}=\begin{bmatrix}
q_{11} & q_{12} & q_{13} & q_{14}\\
q_{21} & q_{22} & q_{23} & q_{24}\\
q_{31} & q_{32} & q_{33} & q_{34}\\
q_{41} & q_{42} & q_{43} & q_{44}
\end{bmatrix} \Longleftrightarrow \begin{bmatrix}
a^{2} & ab & ac & ad\\
ab & b^{2} & bc & bd\\
ac & bc & c^{2} & cd\\
ad & bd & cd & d^{2}
\end{bmatrix}
$$

边坍缩的误差代价

对于edge collapse(contract)操作$$( v_{1} ,\ v_{2}) \ \rightarrow \ v’$$

假设对于某一坍缩边$(v1, v2)$,其收缩后顶点变为 $v’$,我们定义顶点 $v’$ 的误差矩阵为两个合并前的顶点误差矩阵的和:$Q_1+Q_2$,并且可以得到误差代价$\Delta (v’)$
$$\Delta (v’)=v’^{T}( Q_{1} +Q_{2})v’=v’^{T} Q_{e}v’$$
如果我们希望坍缩操作对模型本身的影响最小,那么就是希望 $\Delta (v’)$ 取到最小值。

通过将前面的误差矩阵的矩阵形式替换$\Delta (v’)$,我们可以得到
$$
\Delta (v’)=v’^{T}\begin{bmatrix}
q_{11} & q_{12} & q_{13} & q_{14}\\
q_{21} & q_{22} & q_{23} & q_{24}\\
q_{31} & q_{32} & q_{33} & q_{34}\\
q_{41} & q_{42} & q_{43} & q_{44}
\end{bmatrix}v’
$$
$\Delta (v’)$写成向量$v’( x,\ y,\ z)$并用齐次坐标形式做矩阵乘法化简后可得到
$$
\begin{array}{l}
\Delta (v’)=\begin{bmatrix}
x & y & z & 1
\end{bmatrix}\begin{bmatrix}
q_{11} & q_{12} & q_{13} & q_{14}\\
q_{21} & q_{22} & q_{23} & q_{24}\\
q_{31} & q_{32} & q_{33} & q_{34}\\
q_{41} & q_{42} & q_{43} & q_{44}
\end{bmatrix}\begin{bmatrix}
x\\
y\\
z\\
1
\end{bmatrix}\\
=\begin{bmatrix}
q_{11} x+q_{21} y+q_{31} z+q_{41} & q_{12} x+q_{22} y+q_{32} z+q_{42} & q_{13} x+q_{23} y+q_{33} z+q_{43} & q_{14} +q_{24} +q_{34} +q_{44}
\end{bmatrix}\begin{bmatrix}
x\\
y\\
z\\
1
\end{bmatrix}\\
=\left( q_{11} x^{2} +q_{21} xy+q_{31} xz+q_{41} x\right) +\left( q_{12} xy+q_{22} y^{2} +q_{32} yz+q_{42} y\right)\\
\ \ \ \ +\left( q_{13} xz+q_{23} yz+q_{33} z^{2} +q_{43} z\right) +(q_{14} x+q_{24} y+q_{34} z+q_{44} )\\
=\left( q_{11} x^{2} +q_{22} y^{2} +q_{33} z^{2}\right) +2(q_{14} x+q_{24} y+q_{34} z)+2(q_{12} xy+q_{13} xz+q_{23} yz)+q_{44}
\end{array}
$$
上面的计算过程的最后一步其实是有一个隐式的替换。同时因为误差矩阵Q是一个对称矩阵,所以$q_{ij}=q_{ji}$(例如 $q_{14}=q_{41}$),所以可以合并同类项。
幸运的是,由于$\Delta (v’)$的表达式最终是一个二次项形式,我们可以通过找到一个 $\Delta (v’)’$为0的值来找出它的最小值。即令它对x,y,z的一阶偏导都为0
$$\frac{\partial \Delta (v’)}{\partial x} =\frac{\partial \Delta (v’)}{\partial y} =\frac{\partial \Delta (v’)}{\partial z} =0$$
这里演示$\frac{\partial \Delta (v’)}{\partial x}$的推导
$$
\begin{array}{l}
\frac{\partial \Delta (v’)}{\partial x} =2q_{11} x+2q_{12} y+2q_{13} z+2q_{14}\\
=2\begin{bmatrix}
q_{11} & q_{12} & q_{13} & q_{14}
\end{bmatrix}\begin{bmatrix}
x\\
y\\
z\\
1
\end{bmatrix}\\
\therefore \begin{bmatrix}
q_{11} & q_{12} & q_{13} & q_{14}
\end{bmatrix} \ v’=0
\end{array}
$$
同理可以应用于对y,z的偏微分。故可得
$$
\begin{bmatrix}
q_{11} & q_{12} & q_{13} & q_{14}\\
q_{21} & q_{22} & q_{23} & q_{24}\\
q_{31} & q_{32} & q_{33} & q_{34}\\
0 & 0 & 0 & 1
\end{bmatrix}\begin{bmatrix}
x\\
y\\
z\\
1
\end{bmatrix} =\begin{bmatrix}
0\\
0\\
0\\
1
\end{bmatrix}
$$
因为最下面一行不用处理,将左边的4x4矩阵还可以拆成一个3x3矩阵和3x1矩阵分别乘以 $v’$ 所以式子还可以写成
$$
\begin{array}{l}
\begin{bmatrix}
q_{11} & q_{12} & q_{13}\\
q_{21} & q_{22} & q_{23}\\
q_{31} & q_{32} & q_{33}
\end{bmatrix}\begin{bmatrix}
x\\
y\\
z
\end{bmatrix} +\begin{bmatrix}
q_{14}\\
q_{24}\\
q_{34}
\end{bmatrix}\begin{bmatrix}
1\\
1\\
1
\end{bmatrix} =\begin{bmatrix}
0\\
0\\
0
\end{bmatrix}\\
\begin{bmatrix}
q_{11} & q_{12} & q_{13}\\
q_{21} & q_{22} & q_{23}\\
q_{31} & q_{32} & q_{33}
\end{bmatrix}\begin{bmatrix}
x\\
y\\
z
\end{bmatrix} =-\begin{bmatrix}
q_{14}\\
q_{24}\\
q_{34}
\end{bmatrix}\\
\\
v’=\begin{bmatrix}
q_{11} & q_{12} & q_{13}\\
q_{21} & q_{22} & q_{23}\\
q_{31} & q_{32} & q_{33}
\end{bmatrix}^{-1}\begin{bmatrix}
-q_{14}\\
-q_{24}\\
-q_{34}
\end{bmatrix}
\end{array}
$$

由前面推导的,用平面方程表示的误差矩阵$Q_p$,我们可以得到
$$
v’=\begin{bmatrix}
a^{2} & ab & ac\\
ab & b^{2} & bc\\
ac & bc & c^{2}
\end{bmatrix}^{-1}\begin{bmatrix}
-ad\\
-bd\\
-cd
\end{bmatrix}
$$
这个$v’$也即理论上我们想要找的边坍缩后,最优点的位置。不过这个理论的位置不一定一直存在。例如系数逆矩阵如果不存在逆矩阵,就无法求解。这个时候我们就只能使用备选方案,
在$v_{1} ,\ v_{2}$ 和$( v_{1} +\ v_{2}) /2$三个位置中选择一个最合适的点的位置进行坍缩。
有了最优点,再代回
$$\Delta (v’)=v’^{T}( Q_{1} +Q_{2})v’=v’^{T} Q_{e}v’$$
我们就可以知道某一条边的坍缩代价

坍缩策略

知道了每条边的坍缩代价之后,我们要想尽可能的保留模型原始形状的同时简化边,只要采用贪心的策略不断删去坍缩代价最小的边即可。虽然这个坍缩操作是local operation,不过仍需记得更新相应的相邻点,边,面的坍缩代价。

算法实现细节

基本步骤:
1.首先遍历mesh上所有面,对每一个面计算它的误差矩阵(这与顶点位置无关,所以可以先求)即
$$
Q_p=\begin{bmatrix}
a\\
b\\
c\\
d
\end{bmatrix}\begin{bmatrix}
a & b & c & d
\end{bmatrix}
$$
平面法线$\vec{n}( a,\ b,\ c)$可以由平面法线f->normal()得到。平面的表达式中的d,可以将平面上任意一个顶点的坐标代入平面表达式$ax+by+cy+d=0$中解得。最后在算Q的时候,也不必逐项相乘,如果可以直接使用外积函数自然最为方便

1
f->quadric = outer(Vector4D(n, d), Vector4D(n, d));

2.遍历mesh上所有顶点,对每一个顶点计算它的1-ring相邻面的误差矩阵,直接累加即可。

1
2
3
4
5
6
7
8
9
10
11
for (VertexIter v = mesh.verticesBegin(); v != mesh.verticesEnd(); ++v) {
v->quadric.zero();
HalfedgeCIter h = v->halfedge(); // get one of the outgoing halfedges of the vertex
do {
FaceCIter f = h->face();
// Add up all quadric of vert's 1-ring faces
v->quadric += f->quadric;

h = h->twin()->next(); // move to the next outgoing halfedge of the vertex.
} while (h != v->halfedge()); // keep going until we're back at the beginning
}

3.遍历mesh上所有边,对于每条边,将它的两个端点的误差矩阵和作为边的误差矩阵。然后计算它的最优坍缩点及其代价。
由于我们的坍缩操作是贪心的,每次只取坍缩代价最小的那条边。每次坍缩完成后,因为更新了顶点的位置和周围的1-ring点关联。所以我们会局部的更新相邻的所有涉及到变化的顶点的误差矩阵,再相应的更新面和边的误差矩阵,以便我们下一次迭代能得到最新的坍缩代价。这里推荐使用优先队列来存储每条边的坍缩记录,这样将复杂度降为O(log(N))。而使用一个额外的CollapseRecord结构的原因有两个:1.可以事先将最优点和代价存储在记录中,并持有边的引用。这样优先队列在排序的时候会减少存储空间。2.记录下坍缩的原始顶点可以帮助后面可能进行的还原操作。
这里坍缩的过程也需要小心。在处理每一条需要坍缩的边的时候,要先将其两个端点相关的1-ring领边的坍缩记录从队列中删去(因为后面需要更新),再进行坍缩。坍缩完成之后更新新顶点的1-ring相邻面。以及更新新的顶点的quadric误差矩阵。最后将新的坍缩记录插回优先队列完成更新。

1
2
3
4
5
6
7
8
// Insert any edge touching the new vertex into the queue, creating new edge records for each of them.
h = newV->halfedge();
do {
h->edge()->record = CollapseRecord(h->edge());
queue.insert(h->edge()->record);

h = h->twin()->next();
} while (h != newV->halfedge());

这样才算是完成了一次坍缩。这个过程可能需要多次进行才能将mesh化简实现down sampling,也可以设定一个目标的边的数量进行循环。

示例效果如下

参考

[1] Scotty3D - Simplification
[2] Article Study : Surface Simplification using Quadric Error Metrics
[3] 三维网格精简算法(Quadric Error Metrics)附源码
[4] Edge Collapse
[5] Surface Simplification Using Quadric Error Metrics
[6] Mesh Optimization
[7] Polygonal Mesh - Data Structure and Processing