网格细分的原理其实并不难理解，它的难点主要在于如何实现。在看过无数有原理无代码的博客后，终于决定写一写我的实现方法，并附上代码供大家参考。c++写的可能比较笨拙，望见谅。
1.半边数据结构
很好理解，就是把网格的每一条边分成两个半边，半边是有方向的同一条边的两个半边方向相反。并且一条边是属于两个面，则半边完全属于一个面。
综合上述，得到半边的数据结构为一条半边的起始顶点origin，这条半边指向的下一条半边next，这条半边对应的另外半边opposite，和这条半边所在的面IncFace。
-
typedef
struct HalfEdge//半边结构
-
-
{
-
-
int origin;
-
-
struct HalfEdge* next;
-
-
struct HalfEdge* opposite;
-
-
int IncFace;
-
-
}HalfEdge;
2.网格细分算法Loop subdivision
我只写了三角形网格的细分，即实现了Loop subdivision算法。
首先，我们应该清楚这个算法是怎样细分的：
图中上面两种情况适用于内部点和内部边产生新点的更新，下面两种情况适用于边界点和边界产生新点的更新。
由于我只细分了闭合的三维图形，所以代码只实现了第一种情况，就不讨论第二种了。
细分过程中，
- 每一条边都会产生一个新的顶点，这个顶点由该边所在的两个面中所有顶点决定，即v = 3/8*(v1 + v3) + 1/8*(v0 + v2)
- 每一个原来的顶点都会更新自己的坐标，这个更新由所有与它相邻的原顶点决定，即
，其中
所以，每次细分要做的就是处理每一个旧顶点和每一条边，生成新的顶点，再将新的顶点连接成半边和面。
算法思路：
1.读入图形文件，分别定义数据结构保存图形的点，线，面。定义半边数据结构，将所有边转换为半边。然后开始细分。
2.更新旧顶点：第一步，遍历所有半边，找到一条以该顶点为起始顶点的半边。第二步，通过这条半边找到它的next半边，next半边的origin为v1。再找next半边的next半边 的opposite半边，进入下一个面。第三步，重复第二步直到找到的origin点vi等于v1，则已经找到所有vi，可计算出旧顶点的更新后坐标。
3.在每条旧边上生成新顶点：遍历每一条半边，通过半边和半边的opposite边的next与origin找到两个面的所有顶点，计算出新的顶点。
4.已经得到了所有新的顶点，连接它们得到新的半边和面。画出新图，细分完成。
大概就是这样了，一些地方不知道是否描述清楚，不清楚的请通过代码自行理解。谢谢~
附：（写于ubuntu16.04，g++编译）
-
#include <GL/glut.h>
-
#include <iostream>
-
#include <stdio.h>
-
#include <stdlib.h>
-
#include <math.h>
-
#include <assert.h>
-
#include <string.h>
-
#include<vector>
-
-
using
namespace
std;
-
-
#define PI 3.1415926536
-
-
typedef
struct Vertex
-
{
-
float x,y,z;
-
}Vertex;
-
-
typedef
struct Face
-
{
-
int num;
-
int order[
3];
-
}Face;
-
-
typedef
struct HalfEdge//半边结构
-
{
-
int origin;
-
struct HalfEdge* next;
-
struct HalfEdge* opposite;
-
int IncFace;
-
}HalfEdge;
-
-
typedef
struct Map
-
{
-
int vs,ve,e;
-
-
}Map;
-
-
static
char *filename=
"cube.off";
-
-
vector<Vertex> vertex;
-
vector<Face> face;
-
vector<HalfEdge*> edge;
-
int e_num;
-
int n_node,n_face,n_edge;
-
int width=
800;
-
int height=
800;
-
-
int readoff(const char* filename)
-
{
-
FILE *fp;
-
-
if(!(fp=fopen(filename,
"r")))
-
{
-
fprintf(
stderr,
"Open fail");
-
return
0;
-
}
-
-
char buffer[
1024];
-
if(fgets(buffer,
1023,fp))
-
{
-
if(!
strstr(buffer,
"OFF"))
-
{
-
printf(
"It's not a OFF FILE");
-
return
0;
-
}
-
-
if(fgets(buffer,
1023,fp))
-
{
-
sscanf(buffer,
"%d %d %d",&n_node,&n_face,&n_edge);
-
-
for(
int i=
0;i<n_node;i++)
-
{
-
Vertex ver;
-
fgets(buffer,
1023,fp);
-
sscanf(buffer,
"%f%f%f",&ver.x,&ver.y,&ver.z);
-
vertex.push_back(ver);
-
}
-
for(
int i=
0;i<n_face;i++)
-
{
-
fgets(buffer,
1023,fp);
-
Face f;
-
sscanf(buffer,
"%d%d%d%d",&f.num,&f.order[
0],&f.order[
1],&f.order[
2]);
-
face.push_back(f);
-
}
-
}
-
-
}
-
/*
-
for(int i=0;i<n_node;i++)
-
{
-
cout<<vertex[i].x<<" "<<vertex[i].y<<" "<<vertex[i].z<<endl;
-
}
-
*/
-
}
-
void initEdge()//生成半边存入vector
-
{
-
int
map[n_node][n_node]={
0};
-
for(
int i=
0;i<n_node;i++)
-
{
-
for(
int j=
0;j<n_node;j++)
-
{
-
map[i][j]=
-1;
-
}
-
}
-
cout<<
"map="<<
sizeof(
map[
0])/
sizeof(
int)<<
endl;
-
e_num=
0;
-
for(
int i=
0;i<n_face;i++)
-
{
-
HalfEdge* edge1=
new HalfEdge();
-
HalfEdge* edge2=
new HalfEdge();
-
HalfEdge* edge3=
new HalfEdge();
-
-
edge1->origin=face[i].order[
0];
-
edge2->origin=face[i].order[
1];
-
edge3->origin=face[i].order[
2];
-
-
edge1->next=edge2;
-
edge2->next=edge3;
-
edge3->next=edge1;
-
-
HalfEdge* tmpe=
new HalfEdge();
-
if(
map[face[i].order[
1]][face[i].order[
0]]!=
-1)
-
{
-
tmpe=edge[
map[face[i].order[
1]][face[i].order[
0]]];
-
edge1->opposite=tmpe;
-
tmpe->opposite=edge1;
-
}
-
else
-
{
-
edge1->opposite=
NULL;
-
map[face[i].order[
0]][face[i].order[
1]]=e_num;
-
}
-
e_num++;
-
if(
map[face[i].order[
2]][face[i].order[
1]]!=
-1)
-
{
-
tmpe=edge[
map[face[i].order[
2]][face[i].order[
1]]];
-
edge2->opposite=tmpe;
-
tmpe->opposite=edge2;
-
}
-
else
-
{
-
edge2->opposite=
NULL;
-
map[face[i].order[
1]][face[i].order[
2]]=e_num;
-
}
-
e_num++;
-
if(
map[face[i].order[
0]][face[i].order[
2]]!=
-1)
-
{
-
tmpe=edge[
map[face[i].order[
0]][face[i].order[
2]]];
-
edge3->opposite=tmpe;
-
tmpe->opposite=edge3;
-
}
-
else
-
{
-
edge3->opposite=
NULL;
-
map[face[i].order[
2]][face[i].order[
0]]=e_num;
-
}
-
e_num++;
-
-
edge1->IncFace=i;
-
edge2->IncFace=i;
-
edge3->IncFace=i;
-
-
edge.push_back(edge1);
-
edge.push_back(edge2);
-
edge.push_back(edge3);
-
}
-
n_edge=edge.size();
-
-
//cout<<n_edge<<endl;
-
//for(int i=0;i<n_edge;i++)
-
//cout<<edge[i].origin<<" "<<edge[i].next->origin<<" "<<edge[i].next->next->origin<<" "<<edge[i].IncFace<<endl;
-
//cout<<edge[i]->origin<<" "<<edge[i]->opposite->origin<<endl;
-
}
-
-
HalfEdge* findOriginEdge(int v)//找到从该定点出发的一条半边
-
{
-
for(
int k=
0;k<n_edge;k++)
-
{
-
if(edge[k]->origin==v)
-
return edge[k];
-
}
-
return
NULL;
-
}
-
-
void subdivide()
-
{
-
vector<Vertex> vertex2;
-
vector<Face> face2;
-
vector<HalfEdge*> edge2;
-
HalfEdge* he=
new HalfEdge();
-
int n;
-
float p_sumx,p_sumy,p_sumz;
-
float px,py,pz;
-
float beta;
-
cout<<
"细分开始"<<
endl;
-
for(
int i=
0;i<n_node;i++)
//旧点更新
-
{
-
he=findOriginEdge(i);
-
-
if(he!=
NULL)
-
{
-
n=
0;
-
p_sumx=
0;
-
p_sumy=
0;
-
p_sumz=
0;
-
HalfEdge* e=
new HalfEdge();
-
e=he->next;
-
int p0=e->origin;
-
-
while(e->next->origin!=p0)
-
{
-
n++;
-
p_sumx+=vertex[e->next->origin].x;
-
p_sumy+=vertex[e->next->origin].y;
-
p_sumz+=vertex[e->next->origin].z;
-
HalfEdge* te=
new HalfEdge();
-
te=e->next->opposite;
-
e=te->next;
-
}
-
n++;
-
p_sumx+=vertex[p0].x;
-
p_sumy+=vertex[p0].y;
-
p_sumz+=vertex[p0].z;
-
beta=
1/(
double)n*(
0.625-
pow(
0.375+
0.25*
cos(
2*PI/n),
2));
-
-
px=(
1-n*beta)*vertex[i].x+beta*p_sumx;
-
py=(
1-n*beta)*vertex[i].y+beta*p_sumy;
-
pz=(
1-n*beta)*vertex[i].z+beta*p_sumz;
-
-
Vertex v;
-
v.x=px;
-
v.y=py;
-
v.z=pz;
-
vertex2.push_back(v);
-
}
-
}
-
int map1[n_node][n_node]={
0};
-
-
cout<<
"map1="<<
sizeof(map1[
0])/
sizeof(
int)<<
endl;
-
float qx,qy,qz;
-
-
for(
int i=
0;i<n_edge;i++)
//新点生成
-
{
-
if(!map1[edge[i]->origin][edge[i]->next->origin])
-
{
-
int p=edge[i]->origin;
-
int pi=edge[i]->next->origin;
-
int pi1=edge[i]->next->next->origin;
-
int pi0=edge[i]->opposite->next->next->origin;
-
qx=
0.375*(vertex[p].x+vertex[pi].x)+
0.125*(vertex[pi1].x+vertex[pi0].x);
-
qy=
0.375*(vertex[p].y+vertex[pi].y)+
0.125*(vertex[pi1].y+vertex[pi0].y);
-
qz=
0.375*(vertex[p].z+vertex[pi].z)+
0.125*(vertex[pi1].z+vertex[pi0].z);
-
-
Vertex v;
-
v.x=qx;
-
v.y=qy;
-
v.z=qz;
-
vertex2.push_back(v);
-
-
map1[edge[i]->origin][edge[i]->next->origin]=vertex2.size()
-1;
-
map1[edge[i]->next->origin][edge[i]->origin]=vertex2.size()
-1;
-
}
-
}
-
/*
-
cout<<"新点"<<endl;
-
for(int i=0;i<vertex2.size();i++)
-
{
-
cout<<vertex2[i].x<<" "<<vertex2[i].y<<" "<<vertex2[i].z<<endl;
-
}
-
*/
-
for(
int i=
0;i<n_face;i++)
//新面
-
{
-
int a,b,c,d,e,f;
-
a=face[i].order[
0];
-
b=face[i].order[
1];
-
c=face[i].order[
2];
-
d=map1[a][b];
-
e=map1[b][c];
-
f=map1[a][c];
-
-
Face f2;
-
f2.num=
3;
-
-
f2.order[
0]=a;
-
f2.order[
1]=d;
-
f2.order[
2]=f;
-
face2.push_back(f2);
-
-
f2.order[
0]=d;
-
f2.order[
1]=b;
-
f2.order[
2]=e;
-
face2.push_back(f2);
-
-
f2.order[
0]=d;
-
f2.order[
1]=e;
-
f2.order[
2]=f;
-
face2.push_back(f2);
-
-
f2.order[
0]=f;
-
f2.order[
1]=e;
-
f2.order[
2]=c;
-
face2.push_back(f2);
-
}
-
/*
-
cout<<"新面"<<endl;
-
for(int i=0;i<face2.size();i++)
-
{
-
cout<<face2[i].order[0]<<" "<<face2[i].order[1]<<" "<<face2[i].order[2]<<endl;
-
}
-
*/
-
-
n_face=face2.size();
-
n_node=vertex2.size();
-
cout<<n_node<<
" "<<n_face<<
endl;
-
-
int map2[n_node][n_node]={
0};
-
//int * map2=new int[n_node][n_node];
-
for(
int i=
0;i<n_node;i++)
-
{
-
-
for(
int j=
0;j<n_node;j++)
-
{
-
map2[i][j]=
-1;
-
}
-
}
-
cout<<
"map2="<<
sizeof(map2[
0])/
sizeof(
int)<<
endl;
-
-
e_num=
0;
-
for(
int i=
0;i<n_face;i++)
//新边
-
{
-
HalfEdge* edge4=
new HalfEdge();
-
HalfEdge* edge5=
new HalfEdge();
-
HalfEdge* edge6=
new HalfEdge();
-
-
edge4->origin=face2[i].order[
0];
-
edge5->origin=face2[i].order[
1];
-
edge6->origin=face2[i].order[
2];
-
-
edge4->next=edge5;
-
edge5->next=edge6;
-
edge6->next=edge4;
-
-
HalfEdge* tmpe=
new HalfEdge();
-
if(map2[face2[i].order[
1]][face2[i].order[
0]]!=
-1)
-
{
-
tmpe=edge2[map2[face2[i].order[
1]][face2[i].order[
0]]];
-
edge4->opposite=tmpe;
-
tmpe->opposite=edge4;
-
}
-
else
-
{
-
edge4->opposite=
NULL;
-
map2[face2[i].order[
0]][face2[i].order[
1]]=e_num;
-
}
-
e_num++;
-
if(map2[face2[i].order[
2]][face2[i].order[
1]]!=
-1)
-
{
-
tmpe=edge2[map2[face2[i].order[
2]][face2[i].order[
1]]];
-
edge5->opposite=tmpe;
-
tmpe->opposite=edge5;
-
}
-
else
-
{
-
edge5->opposite=
NULL;
-
map2[face2[i].order[
1]][face2[i].order[
2]]=e_num;
-
}
-
e_num++;
-
if(map2[face2[i].order[
0]][face2[i].order[
2]]!=
-1)
-
{
-
tmpe=edge2[map2[face2[i].order[
0]][face2[i].order[
2]]];
-
edge6->opposite=tmpe;
-
tmpe->opposite=edge6;
-
}
-
else
-
{
-
edge6->opposite=
NULL;
-
map2[face2[i].order[
2]][face2[i].order[
0]]=e_num;
-
}
-
e_num++;
-
-
edge4->IncFace=i;
-
edge5->IncFace=i;
-
edge6->IncFace=i;
-
-
edge2.push_back(edge4);
-
edge2.push_back(edge5);
-
edge2.push_back(edge6);
-
}
-
n_edge=edge2.size();
-
/*
-
cout<<"新边"<<endl;
-
for(int i=0;i<edge2.size();i++)
-
{
-
cout<<edge2[i]->origin<<" "<<edge2[i]->next->origin<<" "<<edge2[i]->IncFace<<endl;
-
}
-
*/
-
vertex.assign(vertex2.begin(),vertex2.end());
-
-
face.assign(face2.begin(),face2.end());
-
-
edge.assign(edge2.begin(),edge2.end());
-
-
cout<<
"完成一次细分"<<
endl;
-
cout<<n_node<<
" "<<n_edge<<
" "<<n_face<<
endl;
-
}
-
-
void display()
-
{
-
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
//清除颜色和深度缓存
-
glMatrixMode(GL_MODELVIEW);
-
glLoadIdentity();
//重置当前的模型观察矩阵
-
glPushMatrix();
-
glTranslatef(
0.0f,
0.0f,
-5.0f);
-
-
glRotatef(
30,
0.0f,
1.0f,
0.0f);
//饶轴旋转
-
glRotatef(
30,
1.0f,
0.0f,
0.0f);
-
glColor3f(
0.5f,
0.5f,
0.5f);
//灰色
-
-
glBegin(GL_TRIANGLES);
-
for(
int i=
0;i<n_face;i++)
-
{
-
glVertex3f(vertex[face[i].order[
0]].x,vertex[face[i].order[
0]].y,vertex[face[i].order[
0]].z);
-
glVertex3f(vertex[face[i].order[
1]].x,vertex[face[i].order[
1]].y,vertex[face[i].order[
1]].z);
-
glVertex3f(vertex[face[i].order[
2]].x,vertex[face[i].order[
2]].y,vertex[face[i].order[
2]].z);
-
}
-
glEnd();
-
glPopMatrix();
-
glutSwapBuffers();
-
}
-
-
void keyboard(unsigned char key, int x, int y)
-
{
-
switch(key)
-
{
-
case
'1':
-
glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
-
break;
-
case
'2':
-
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
-
break;
-
case
'3':
-
glPolygonMode(GL_FRONT_AND_BACK, GL_POINT);
-
break;
-
case
'w':
-
subdivide();
-
break;
-
}
-
glutPostRedisplay();
-
}
-
-
void reshape(int w, int h) {
-
//定义视口大小
-
glViewport(
0,
0, (GLsizei) w, (GLsizei) h);
-
//投影显示
-
glMatrixMode(GL_PROJECTION);
-
//坐标原点在屏幕中心
-
glLoadIdentity();
-
//操作模型视景
-
gluPerspective(
60.0, (GLfloat) w/(GLfloat) h,
1.0,
20.0);
-
glMatrixMode(GL_MODELVIEW);
-
}
-
-
int main(int argc,char **argv)
-
{
-
readoff(filename);
-
initEdge();
-
-
glutInit(&argc,argv);
-
glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB|GLUT_DEPTH);
-
glutInitWindowSize(width,height);
-
glutInitWindowPosition(
100,
100);
-
glutCreateWindow(
"loop");
-
glutReshapeFunc(reshape);
-
glutDisplayFunc(display);
-
glutIdleFunc(display);
-
glutKeyboardFunc(keyboard);
-
glutMainLoop();
-
return
0;
-
}
图形文件cube.off
-
OFF
-
8 12 0
-
-0
.5
-0
.5
-0
.5
-
0
.5
-0
.5
-0
.5
-
-0
.5 0
.5
-0
.5
-
0
.5 0
.5
-0
.5
-
-0
.5
-0
.5 0
.5
-
0
.5
-0
.5 0
.5
-
-0
.5 0
.5 0
.5
-
0
.5 0
.5 0
.5
-
3 4 5 6
-
3 6 5 7
-
3 5 1 7
-
3 7 1 3
-
3 0 3 1
-
3 0 2 3
-
3 4 6 2
-
3 4 2 0
-
3 6 7 2
-
3 3 2 7
-
3 5 4 0
-
3 5 0 1
2459
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
