【几何处理】基于二次误差测量的网格简化（二）

最新推荐文章于 2021-04-01 17:17:27 发布

旅程TSH

最新推荐文章于 2021-04-01 17:17:27 发布

阅读量2k

点赞数 4

分类专栏：数字几何处理

本文链接：https://blog.csdn.net/u010669231/article/details/104450059

版权

数字几何处理专栏收录该内容

5 篇文章 11 订阅

订阅专栏

接上文，Hope的边收缩操作可推广为一般的顶点合并变换来描述 $(v_{1},v_{2})\rightarrow v$ ,其含义意是将两个顶点 $v_{1}v_{2}$ 移到一组新的位置v，将连向 $v_{1}v_{2}$ 的所有边都连向v,并删除所有退化的边和面片。

基本的思想就是不断地收缩边，直到达到我们所需要的三角形数量，我们收缩的边越多，网格越简单。那么，我们应该收缩哪些边，当我们收缩边时，我们应该在哪里放置新顶点？Garland和Heckbert提出了一种简单的贪心方案，该方案在实践中效果很好，并且是当今许多网格简化工具的基础。

Garland和Heckbert进入了二次误差度量来刻画每一个顶点移动后的误差，对表面上每一个顶点 $v_{a}$ 均有许多三角面片与之相邻，记 $plane(v_{a})$ 为这些三角形所在平面方程所构成的集合，即

$plane(v_{a})={(a,b,c,d)|ax+by+cz+d=0(a^{2}+b^{2}+c^{2}=1)}$

则我们采用如下的二次函数来度量va移动到v产生的误差

$\Delta (v_{a}\rightarrow v)=\sum_{plane(v_{a})}(pv^{T})^{2}$

其中 $v=(x,y,z,1)$ 为齐次坐标。展开上式得到

$\Delta (v_{a}\rightarrow v)=\sum_{plane(v_{a})}(pv^{T})^{2}=v(\sum_{plane(v_{a})}K_{p})v^{T}=vQ(v_{a})v^{T}$ 。

我们知道pv=0就是平面方程，如果这个点没有变化，我们把v带进该误差公式，得到的误差就等于0。由于我们在考虑误差时不考虑是在平面上还是平面下，我们我们对该公式进行平方。将点va周围一圈三角形的误差加起来，

上式中

$K_{p}=p^{T}p=\begin{bmatrix} a^2 &ab &ac &ad \\ ab &b^2 &bc &bd \\ ac&bc &c^{2} &cd \\ ad&bd &cd &d^{2} \end{bmatrix}$

	for (auto f = mesh.facesBegin(); f != mesh.facesEnd(); f++)
	{
		double d = -dot(f->normal(), f->halfedge()->vertex()->position);
		Vector4D v = Vector4D(f->normal(), d);
		f->quadric = outer(v, v);
	}

$Q(v_{a})=\sum_{plane(v_{a})}K_{p}$

	for (auto v = mesh.verticesBegin(); v != mesh.verticesEnd(); v++)
	{
		auto adjFs = v->AdjFaces();
		v->quadric.zero();
		for (auto f : adjFs)
			v->quadric += f->quadric;
	}

这样，对每一组顶点v_{a},在预处理时，我们均可按上述方法计算矩阵 $Q(v_{a})$ ,进而就可对其移动进行误差度量了。但由于每次合并时，须同时移动两点，故必须考虑同时移动多个点后产生的误差。

Garland和Heckbert简单地采用加法规则来刻画多点移动形成的误差，对点对合并 $(v_{1},v_{2})\rightarrow v$ ,其误差为 $\Delta (v)=\Delta (v_{1}\rightarrow v)+\Delta (v_{2}\rightarrow v)=v(Q(v_{1})+Q(v_{2}))v^{T}=vQv^{T}$ 。

Matrix4x4 quadricE = edge->halfedge()->vertex()->quadric + edge->halfedge()->twin()->vertex()->quadric;

计算误差代码：

Vector4D u = Vector4D(optimalPoint, 1);
score = dot((quadricE*u), u);

如果收缩边的话，我们应该将新点放在哪里，才能使得误差最小。换句话说，应该放哪里才能使 $x^{T}Kx$ 最小化？其实就是求解 $Ku=0$ 。对于用齐次坐标表示的最优位置u，我们可以根据齐次坐标的w分量为1来简化这个过程。在经过一些处理后，我们将这个方程简化为一个更小的3*3的线性方程 $Ax=b$ ,这里

$A=\begin{bmatrix} a^2 &ab &ac \\ ab &b^2 &bc \\ ac&bc &c^{2}\end{bmatrix}$ $b=\begin{bmatrix} -ad \\ -bd \\ -cd\end{bmatrix}$ 。

Matrix3x3 A;
for (size_t x = 0; x <= 2; x++)
{
	for (size_t y = 0; y <= 2; y++)
	{
		A(x, y) = quadricE(x, y);
	}
}
Vector3D w(quadricE(0, 3), quadricE(1, 3), quadricE(2, 3));
Vector3D b = -w;

如果A不可逆的话我们就简单的将v设为 $\frac{v_{1}+v_{2}}{2}$ ,如果有解通过 $x=A^{-1}b$ 求解即可

	if (abs(A.det()) < 0.000001) {
		optimalPoint = edge->centroid();
	}
	else
	{
		optimalPoint = A.inv()*b;
	}

为了使算法的误差最小，我们将所有边放入一个优先级队列，该队列是根据上文计算出的误差score来排序的，每次进行边收缩的时候，选取队首元素即可。

	MutablePriorityQueue<EdgeRecord> queue;
	for (auto e = mesh.edgesBegin(); e != mesh.edgesEnd(); e++)
	{
		e->record = EdgeRecord(e);
		queue.insert(e->record);
	}

下面是EdgeRecord的定义

class EdgeRecord {
 public:
  EdgeRecord() {}
  EdgeRecord(EdgeIter& _edge);

  EdgeIter edge;
  Vector3D optimalPoint;
  double score;
};

inline bool operator<(const EdgeRecord& r1, const EdgeRecord& r2) {
  if (r1.score != r2.score) {
    return (r1.score < r2.score);
  }

  EdgeIter e1 = r1.edge;
  EdgeIter e2 = r2.edge;
  return &*e1 < &*e2;
}

最后一步概述如下：

1.从队列中获取误差最小的边

2.删除队列中的该边

3.通过加和两端点的误差来计算新误差

4.移除所有该边的邻接边

5.进行边收缩

6.将新误差赋值给新的顶点

7.将新顶点的邻接边全部放入队列，并为每个边设置EdgeRecord。

	size_t targetNum = mesh.nFaces() / 4;
	while (mesh.nFaces()>targetNum)
	{
		EdgeRecord eR = queue.top();
		queue.pop();

		{//remove adjEs' record in queue
			auto adjEs = eR.edge->AdjEdges();
			for (auto adjE : adjEs)
				queue.remove(adjE->record);
		}

		VertexIter newV = mesh.collapseEdge(eR.edge);
		if (!mesh.IsValid(newV, "downsample : collapse an edge fail"))
			return;

		newV->position = eR.optimalPoint;

		{// set adjFs' and newV's quadric
			newV->quadric.zero();
			auto adjFs = newV->AdjFaces();
			for (auto adjF : adjFs) {
				double d = -dot(adjF->normal(), adjF->halfedge()->vertex()->position);
				Vector4D v = Vector4D(adjF->normal(), d);
				adjF->quadric = outer(v, v);
				newV->quadric += adjF->quadric;
			}
		}

		{// set adjVs' quadric
			auto adjVs = newV->AdjVertices();
			for (auto adjV : adjVs) {
				adjV->quadric.zero();
				auto adjFs = adjV->AdjFaces();
				for (auto adjF : adjFs)
					adjV->quadric += adjF->quadric;
			}
		}

		{// set adjEs' record
			auto adjEs = newV->AdjEdges();
			for (auto adjE : adjEs) {
				adjE->record = EdgeRecord(adjE);
				queue.insert(adjE->record);
			}
		}
	}

运行结果：该程序使用的是CMU462课程的Scotty3D程序

未简化：

简化到原面数的1/1.5：

简化到原面数的1/2：

简化到原面数的1/3：

简化到原面数的1/4：

相关程序代码：

EdgeRecord::EdgeRecord(EdgeIter& _edge) : edge(_edge) {
  // TODO: (meshEdit)
  // Compute the combined quadric from the edge endpoints.
  // -> Build the 3x3 linear system whose solution minimizes the quadric error
  //    associated with these two endpoints.
	Matrix4x4 quadricE = edge->halfedge()->vertex()->quadric + edge->halfedge()->twin()->vertex()->quadric;
	Matrix3x3 A;
	for (size_t x = 0; x <= 2; x++)
	{
		for (size_t y = 0; y <= 2; y++)
		{
			A(x, y) = quadricE(x, y);
		}
	}
	Vector3D w(quadricE(0, 3), quadricE(1, 3), quadricE(2, 3));
	Vector3D b = -w;
  // -> Use this system to solve for the optimal position, and store it in
  //    EdgeRecord::optimalPoint.
	if (abs(A.det()) < 0.000001) {
		optimalPoint = edge->centroid();
	}
	else
	{
		optimalPoint = A.inv()*b;
	}
  // -> Also store the cost associated with collapsing this edg in
  //    EdgeRecord::Cost.
	Vector4D u = Vector4D(optimalPoint, 1);
	score = dot((quadricE*u), u);
}

void MeshResampler::downsample(HalfedgeMesh& mesh) {
  // TODO: (meshEdit)
  // Compute initial quadrics for each face by simply writing the plane equation
  // for the face in homogeneous coordinates. These quadrics should be stored
  // in Face::quadric
  // -> Compute an initial quadric for each vertex as the sum of the quadrics
  //    associated with the incident faces, storing it in Vertex::quadric
  // -> Build a priority queue of edges according to their quadric error cost,
  //    i.e., by building an EdgeRecord for each edge and sticking it in the
  //    queue.
  // -> Until we reach the target edge budget, collapse the best edge. Remember
  //    to remove from the queue any edge that touches the collapsing edge
  //    BEFORE it gets collapsed, and add back into the queue any edge touching
  //    the collapsed vertex AFTER it's been collapsed. Also remember to assign
  //    a quadric to the collapsed vertex, and to pop the collapsed edge off the
  //    top of the queue.
	for (auto f = mesh.facesBegin(); f != mesh.facesEnd(); f++)
	{
		double d = -dot(f->normal(), f->halfedge()->vertex()->position);
		Vector4D v = Vector4D(f->normal(), d);
		f->quadric = outer(v, v);
	}

	for (auto v = mesh.verticesBegin(); v != mesh.verticesEnd(); v++)
	{
		auto adjFs = v->AdjFaces();
		v->quadric.zero();
		for (auto f : adjFs)
			v->quadric += f->quadric;
	}

	MutablePriorityQueue<EdgeRecord> queue;
	for (auto e = mesh.edgesBegin(); e != mesh.edgesEnd(); e++)
	{
		e->record = EdgeRecord(e);
		queue.insert(e->record);
	}

	size_t targetNum = mesh.nFaces() / 4;
	while (mesh.nFaces()>targetNum)
	{
		EdgeRecord eR = queue.top();
		queue.pop();

		{//remove adjEs' record in queue
			auto adjEs = eR.edge->AdjEdges();
			for (auto adjE : adjEs)
				queue.remove(adjE->record);
		}

		VertexIter newV = mesh.collapseEdge(eR.edge);
		if (!mesh.IsValid(newV, "downsample : collapse an edge fail"))
			return;

		newV->position = eR.optimalPoint;

		{// set adjFs' and newV's quadric
			newV->quadric.zero();
			auto adjFs = newV->AdjFaces();
			for (auto adjF : adjFs) {
				double d = -dot(adjF->normal(), adjF->halfedge()->vertex()->position);
				Vector4D v = Vector4D(adjF->normal(), d);
				adjF->quadric = outer(v, v);
				newV->quadric += adjF->quadric;
			}
		}

		{// set adjVs' quadric
			auto adjVs = newV->AdjVertices();
			for (auto adjV : adjVs) {
				adjV->quadric.zero();
				auto adjFs = adjV->AdjFaces();
				for (auto adjF : adjFs)
					adjV->quadric += adjF->quadric;
			}
		}

		{// set adjEs' record
			auto adjEs = newV->AdjEdges();
			for (auto adjE : adjEs) {
				adjE->record = EdgeRecord(adjE);
				queue.insert(adjE->record);
			}
		}
	}
}

旅程TSH

关注

4
点赞
踩
9

收藏

觉得还不错? 一键收藏
1
评论
【几何处理】基于二次误差测量的网格简化（二）

接上文，Hope的边收缩操作可推广为一般的顶点合并变换来描述,其含义意是将两个顶点移到一组新的位置v，将连向的所有边都连向v,并删除所有退化的边和面片。基本的思想就是不断地收缩边，直到达到我们所需要的三角形数量，我们收缩的边越多，网格越简单。那么，我们应该收缩哪些边，当我们收缩边时，我们应该在哪里放置新顶点？Garland和Heckbert提出了一种简单的贪心方案，该方案在实践中效...
复制链接

扫一扫