半邊資料結構與網格細分演算法Loop subdivision（附程式碼）
阿新 • • 發佈：2018-11-07
網格細分的原理其實並不難理解，它的難點主要在於如何實現。在看過無數有原理無程式碼的部落格後，終於決定寫一寫我的實現方法，並附上程式碼供大家參考。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.在每條舊邊上生成新頂點
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