CodeRead - Fast quadric mesh simplification
CodeRead - Fast quadric mesh simplification
代码路径见：sp4cerat/Fast-Quadric-Mesh-Simplification: Mesh triangle reduction using quadrics (github.com)
核心思想是基于Qudaric，和[[openmesh-src-decimation]]中介绍的Decimate Algorithm相似。
基本数据结构介绍
// Triangle: 是三角形的定义
// 其中v[3]是三角形三个顶点在std::vector<Vertex> vertices中的索引
// err[0-3]：表示三角形的边0-1,1-2,2-3,每条边利用calculate_error计算得到的误差
// err[4]: 表示上面三个error的最小error，作为该三角形的error
// deleted: 表示该三角形已经被删除；
// dirty：表示在处理的过程中，该三角形有相关的顶点已经参与了collapse；
// n: 法向量
// uvs, material: 目前不关心纹理，材质相关信息；
struct Triangle { int v[3];double err[4];int deleted,dirty,attr;vec3f n;vec3f uvs[3];int material; };
// Vertex：顶点数据结构
// p: 顶点的位置
// tstart: Ref中对应的关联三角形信息的起始位置
// tcount: 该顶点关联的三角形个数
// q: 该顶点对应的Quadric Error Matrix
// border: 判断该顶点是不是位于边界
struct Vertex { vec3f p;int tstart,tcount;SymetricMatrix q;int border;};
// 用来存储每个顶点关联的三角形信息
// tid: 三角形的id，用于从triangles中快速获取对应的三角形
// tvertex: 记录该顶点对应的是三角形中的哪个顶点，0,1,2
struct Ref { int tid,tvertex; };
// 存储所有的三角形
std::vector<Triangle> triangles;
// 存储所有的顶点
std::vector<Vertex> vertices;
// 用来存储每个顶点关联的三角形index
std::vector<Ref> refs;
Simplify_mesh 函数简介
代码的入口函数为simplify_mesh
void simplify_mesh(int target_count, double agreesiveness=7, bool verbose=false)
{
// 1. 标注所有三角形的deleted tag为0
...;
// 2. 迭代loop
for (int iteration = 0; iteration < 100; iteration++)
{
// 2.1 减少的三角形数量满足要求，则退出迭代
if (triangle_count-deleted_triangles <= target_count) break;
// 2.2 每个一定时间更新网格
// update_mesh处理的事情是，将三角形的存储更加紧凑，第一次调用的时候
// 会初始化 Quadrics by Plane & Edge Errors
if (iteration%5 == 0) update_mesh(iteration);
// 2.3 清理三角形的dirty flag
...
// 2.4 神秘公式，用于设定阈值
double threshold = 0.000000001*pow(double(iteration+3),agressiveness);
// 2.5 删除顶点，并且标记删除的三角形
for (int i=0; i<triangles.size(); i++)
{
// 2.5.1 略过部分场景
Triangle &t = triangles[i];
if (t.err[3] > threshold) continue;
if (t.deleted) continue;
if (t.dirty) continue;
// 2.5.2 对三角形的顶点进行遍历
for (int j=0; j<3; ++j)
{
if (t.err[j] < threshold)
{
// 略过一个点在边界上，一个点不在边界上的情况
int i0=t.v[j]; Vertex &v0 = vertices[i0];
int i1=t.v[(j+1)%3]; Vertex &v1 = vertices[i1];
if (v0.border != v1.border) continue;
// 根据线性方程求解，计算塌陷后的顶点的位置
vec3f p;
calculate_error(i0, i1, p);
// 判断塌陷后是不是会出现三角形重叠的情况
deleted0.resize(v0.tcount); // 用来存储v0关联的三角形的法向量
deleted1.resize(v1.tcount); // 用来存储v1关联的三角形的法向量
if (flipped(p,i0,i1,v0,v1,deleted0)) continue;
if (flipped(p,i1,i0,v1,v0,deleted1)) continue;
// 塌陷没有折叠边
v0.p = p;
v0.q = v1.q + v0.q;
int tstart=refs.size();
update_triangles(i0,v0,deleted0,deleted_triangles);
update_triangles(i0,v1,deleted1,deleted_triangles);
// 节约内存的方法
if(tcount<=v0.tcount)
{
// save ram
if(tcount) memcpy(&refs[v0.tstart],&refs[tstart],tcount*sizeof(Ref));
}
else
// append
v0.tstart=tstart;
v0.tcount=tcount;
break;
}
if (triangle_count-deleted_triangles<=target_count)break;
}
}
}
// clean up mesh
compact_mesh();
}
上面主体流程为：
- 初始化计算每个顶点的Quadric Error Metric；
- 初始化计算每条边塌陷后的误差E；
- 遍历所有三角形的边，判断是否满足调整阈值；
- 如果否，计算边塌陷后顶点位置，判断边塌陷后是否会产生折叠；
- 如果否，更新三角形的连接关系，以及边的误差E；
- 如果否，计算边塌陷后顶点位置，判断边塌陷后是否会产生折叠；
其中，没有展开来介绍的函数有：
- update_mesh
- calculate_error
- flipped
- update_triangles
- compact_mesh
Update_mesh函数简介（内带calcuate_error介绍）
关于update_mesh重点需要介绍的功能是初始化Quadrics by plane & edge errors。
void update_mesh(int iteration)
{
.....;
if( iteration == 0 )
{
// Identify boundary : vertices[].border=0,1
std::vector<int> vcount,vids;
loopi(0,vertices.size())
vertices[i].border=0;
// 用来实现判断一个点是不是边界点的逻辑，
// TODO: 此处关于边界点的判断是有问题的。
// 以正六边形为例，每个顶点和正六边形的中心相连,
// 删除其中的一条边，那么中心的点应该是边界点，但是下面的判断，
// 认为中间的点并不是边界点
loopi(0,vertices.size()) // 遍历所有顶点
{
Vertex &v=vertices[i];
vcount.clear();
vids.clear();
loopj(0,v.tcount) // 遍历顶点v周围邻接三角形
{
int k=refs[v.tstart+j].tid;
Triangle &t=triangles[k];
loopk(0,3) // 遍历该三角形的三个顶点
{
int ofs=0,id=t.v[k];
while(ofs<vcount.size())
{
if(vids[ofs]==id)break;
ofs++;
}
if(ofs==vcount.size())
{
vcount.push_back(1);
vids.push_back(id);
}
else
vcount[ofs]++;
}
}
loopj(0,vcount.size()) if(vcount[j]==1)
vertices[vids[j]].border=1;
}
//initialize errors
loopi(0,vertices.size())
vertices[i].q=SymetricMatrix(0.0);
// 计算每个顶点Quadric error matrix
loopi(0,triangles.size())
{
Triangle &t=triangles[i];
vec3f n,p[3];
loopj(0,3) p[j]=vertices[t.v[j]].p;
n.cross(p[1]-p[0],p[2]-p[0]);
n.normalize(); // 计算三角形的法向量
t.n=n;
// 初始化计算每个顶点的Quadric error matrix
loopj(0,3) vertices[t.v[j]].q =
vertices[t.v[j]].q+SymetricMatrix(n.x,n.y,n.z,-n.dot(p[0]));
}
// 计算最小的error，用于和阈值进行比较
loopi(0,triangles.size())
{
// Calc Edge Error
Triangle &t=triangles[i];vec3f p;
loopj(0,3) t.err[j]=calculate_error(t.v[j],t.v[(j+1)%3],p);
t.err[3]=min(t.err[0],min(t.err[1],t.err[2]));
}
}
}
Quadric Error Metric介绍
在Polygon Mesh Processing一书中关于decimator有这么一段内容：
根据对公式7.2求解，可以得到具体的坐标位置x。
SymetricMatrix介绍
代码实现中为了表示Quadric Error Metrix矩阵，引入了SymetricMatrix这个类。
class SymetricMatrix {
// a,b,c 构成平面的法向量，d为平面上的点在法向量上投影距离的负值
SymetricMatrix(double a,double b,double c,double d)
{
m[0] = a*a; m[1] = a*b; m[2] = a*c; m[3] = a*d;
m[4] = b*b; m[5] = b*c; m[6] = b*d;
m[7 ] =c*c; m[8 ] = c*d;
m[9 ] = d*d;
}
};
利用SymetricMatrix计算vertex_error
具体公式如7.1所示。该公式展开后对应的代码实现如下：
double vertex_error(SymetricMatrix q, double x, double y, double z)
{
return q[0]*x*x + 2*q[1]*x*y + 2*q[2]*x*z + 2*q[3]*x + q[4]*y*y
+ 2*q[5]*y*z + 2*q[6]*y + q[7]*z*z + 2*q[8]*z + q[9];
}
利用SymetricMatrix计算vertex_position
该计算公式如7.2所示，该公式求解方式如下：
如果对应的行列式D为0，那么从边的两个端点和中心点处找到vertex error最小的位置。对应的代码如下：
double calculate_error(int id_v1, int id_v2, vec3f &p_result)
{
// compute interpolated vertex
SymetricMatrix q = vertices[id_v1].q + vertices[id_v2].q;
bool border = vertices[id_v1].border & vertices[id_v2].border;
double error=0;
double det = q.det(0, 1, 2, 1, 4, 5, 2, 5, 7);
if ( det != 0 && !border )
{
// q_delta is invertible
p_result.x = -1/det*(q.det(1, 2, 3, 4, 5, 6, 5, 7 , 8)); // vx = A41/det(q_delta)
p_result.y = 1/det*(q.det(0, 2, 3, 1, 5, 6, 2, 7 , 8)); // vy = A42/det(q_delta)
p_result.z = -1/det*(q.det(0, 1, 3, 1, 4, 6, 2, 5, 8)); // vz = A43/det(q_delta)
error = vertex_error(q, p_result.x, p_result.y, p_result.z);
}
else
{
// det = 0 -> try to find best result
vec3f p1=vertices[id_v1].p;
vec3f p2=vertices[id_v2].p;
vec3f p3=(p1+p2)/2;
double error1 = vertex_error(q, p1.x,p1.y,p1.z);
double error2 = vertex_error(q, p2.x,p2.y,p2.z);
double error3 = vertex_error(q, p3.x,p3.y,p3.z);
error = min(error1, min(error2, error3));
if (error1 == error) p_result=p1;
if (error2 == error) p_result=p2;
if (error3 == error) p_result=p3;
}
return error;
}
flipped函数介绍
该函数用来判断，如果一条边被删去之后，这个三角形是否会发生折叠。代码如下：
/// \param[in] p 计算出来新的点的位置
/// \param[in] i0, v0 边的一个点
/// \param[in] i1, v1 边的另一个点
/// \param[in] deleted i0顶点关联的周边三角形的集合，用于存储法向量
bool flipped(vec3f p,int i0,int i1,Vertex &v0,Vertex &v1,std::vector<int> &deleted)
{
loopk(0,v0.tcount)
{
Triangle &t=triangles[refs[v0.tstart+k].tid];
if(t.deleted)continue;
int s=refs[v0.tstart+k].tvertex;
int id1=t.v[(s+1)%3];
int id2=t.v[(s+2)%3];
// 如果三角形的一条边正好是塌陷边，那么这个三角形是即将被删除的三角形
if(id1==i1 || id2==i1) // delete ?
{
deleted[k]=1;
continue;
}
vec3f d1 = vertices[id1].p-p; d1.normalize();
vec3f d2 = vertices[id2].p-p; d2.normalize();
// 如果即将形成的新的三角形的两个边基本同向，返回true
if(fabs(d1.dot(d2))>0.999) return true;
vec3f n;
n.cross(d1,d2);
n.normalize();
deleted[k]=0;
// 如果即将得到的新的三角形的法向量和之前的法向量之间的夹角过大，返回true
if(n.dot(t.n)<0.2) return true;
}
return false;
}
update_triangles函数介绍
// Update triangle connections and edge error after a edge is collapsed
// 更新三角形的连接关系，以及边的error
void update_triangles(int i0,Vertex &v,std::vector<int> &deleted,int &deleted_triangles)
{
vec3f p;
loopk(0,v.tcount)
{
Ref &r=refs[v.tstart+k];
Triangle &t=triangles[r.tid];
// 如果删除，则continue
if(t.deleted)continue;
if(deleted[k])
{
t.deleted=1;
deleted_triangles++;
continue;
}
// 如果没有被删除，重新计算三角形的err，并更新dirty状态
t.v[r.tvertex]=i0;
t.dirty=1;
t.err[0]=calculate_error(t.v[0],t.v[1],p);
t.err[1]=calculate_error(t.v[1],t.v[2],p);
t.err[2]=calculate_error(t.v[2],t.v[0],p);
t.err[3]=min(t.err[0],min(t.err[1],t.err[2]));
refs.push_back(r);
}
}
版权说明
作者： grassofsky
出处： http://www.cnblogs.com/grass-and-moon
本文版权归作者，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出, 原文链接 如有问题， 可邮件（grass-of-sky@163.com）咨询.
