前些天，我一個交大研究生同學給了我一個大作業，讓我嘗試下，內容大致如下：
能夠讀取STL模型文件，獲得其軸對齊包圍盒信息，按照Z軸方向，實用一個無限大平面將該STL模型均勻地剖分5次；5次切分獲得的輪廓線上點存入一個數據庫；軟件能夠從數據庫讀取這5個輪廓，在界面中顯示出來。
這個問題其實還是很簡單的，至於STL文件的解析和讀取在上篇的博文中有所提及，至於剖分取點的內容，其實說白了也就是初中的相似三角形原理，我們給定一個Z為法向的平面，讓它與一個實體求交，在STL模型中，這個實體實際上就是一系列的三角面片，那么我們實際上就是求這個平面和一系列三角面片的交點，很顯然，它們的交集是一條直線（在一般情況下），而且並不是每個三角面片和這個平面都有交線。只有三角形的三個點的Z坐標中有比這個Z值大且有比這個Z值小的時候，才會和這個平面存在交線，很明顯這個時候，這個Z法向的平面把三角面片給切到了。兩點確定一直線，我們只需要求得這個線上的兩個點就行了，這兩個點我們就在三角形的兩邊選取，例如三角形三個點的Z坐標:Z1，Z2，Z2,不妨設Z1<Z，Z2>Z，Z3>Z，那么我們就用點1和點2的直線以及點1和點3的直線與平面求交，這個求交其實也很簡單，即相似三角形的內容，點1和點2組成的直線與z=Z的交點為：x=X1+(Z-Z1)*(X2-X1)/(Z2-Z1)依次類推，可得這兩個點的坐標，那么在所有的三角面片中求得這些點，也就構成了輪廓點，當然，這只是實體布爾運算最簡單的部分，個人拙見，望感興趣的觀眾提出意見。
以下是代碼以及實例，代碼的風格說實話我自己都有點看不下去，也是倉促有靈感寫出來的，包圍盒的算法也很簡單，無非是在這些三角面片中取出最小和最大值，將其圍起來即可
#include "windows.h" #include <GLTools.h> // OpenGL toolkit #include <GLMatrixStack.h> #include <GLFrame.h> #include <GLFrustum.h> #include <GLGeometryTransform.h> #include <iostream> #include <fstream> #include <math.h> #ifdef __APPLE__ #include <glut/glut.h> #else #define FREEGLUT_STATIC #include <GL/glut.h> #endif using namespace std; GLFrame viewFrame; GLFrustum viewFrustum; GLBatch triangleBatch; GLBatch boundaryBatch; GLBatch polyline[100]; GLMatrixStack modelViewMatix; GLMatrixStack projectionMatrix; GLGeometryTransform transformPipeline; GLShaderManager shaderManager; int num; int n,count=0; float* verts; float* vnorms; int tip=1; bool l=true; GLfloat xverts[108]; GLfloat xnormal[108]= { 0,0,-1, 0,0,-1, 0,0,-1, 0,0,-1, 0,0,-1, 0,0,-1, 0,-1,0, 0,-1,0, 0,-1,0, 0,-1,0, 0,-1,0, 0,-1,0, 0,0,1, 0,0,1, 0,0,1, 0,0,1, 0,0,1, 0,0,1, -1,0,0, -1,0,0, -1,0,0, -1,0,0, -1,0,0, -1,0,0, 0,1,0, 0,1,0, 0,1,0, 0,1,0, 0,1,0, 0,1,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0, 1,0,0 }; void RenderScene(); void getpline(float z,int num,GLBatch& batch,GLfloat* verts) { float a,b,c; GLfloat *mp=new GLfloat[6*num]; int p=0,num2=0; float rate; ofstream outfile("1-2.txt", std::ios::app); //outfile.seekp(0,ios::end); outfile<<"\n"; for (int i=0;i<num;i++) { a=verts[i*9+2]; b=verts[i*9+5]; c=verts[i*9+8]; if ((a-z)*(b-z)>0&&(a-z)*(c-z)<0) { rate=(c-z)/(c-a); mp[p]=verts[i*9+6]+rate*(verts[i*9]-verts[i*9+6]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+7]+rate*(verts[i*9+1]-verts[i*9+7]); outfile<<mp[p]<<" "; p++; mp[p]=z; outfile<<mp[p]<<"\n"; p++; rate=(c-z)/(c-b); mp[p]=verts[i*9+6]+rate*(verts[i*9+3]-verts[i*9+6]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+7]+rate*(verts[i*9+4]-verts[i*9+7]); outfile<<mp[p]<<" "; p++; mp[p]=z; outfile<<mp[p]<<"\n"; p++; num2+=2; continue; } if ((a-z)*(c-z)>0&&(a-z)*(b-z)<0) { rate=(b-z)/(b-a); mp[p]=verts[i*9+3]+rate*(verts[i*9]-verts[i*9+3]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+4]+rate*(verts[i*9+1]-verts[i*9+4]); outfile<<mp[p]<<" "; p++; mp[p]=z; outfile<<mp[p]<<"\n"; p++; rate=(b-z)/(b-c); mp[p]=verts[i*9+3]+rate*(verts[i*9+6]-verts[i*9+3]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+4]+rate*(verts[i*9+7]-verts[i*9+4]); outfile<<mp[p]<<" "; p++; mp[p]=z; outfile<<mp[p]<<"\n"; p++; num2+=2; continue; } if ((c>z&&b>z&&a<z)||(c<z&&b<z&&a>z)) { rate=(a-z)/(a-b); mp[p]=verts[i*9]+rate*(verts[i*9+3]-verts[i*9]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+1]+rate*(verts[i*9+4]-verts[i*9+1]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+2]+rate*(verts[i*9+5]-verts[i*9+2]); /*mp[p]=z;*/ outfile<<mp[p]<<"\n"; p++; rate=(a-z)/(a-c); mp[p]=verts[i*9]+rate*(verts[i*9+6]-verts[i*9]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+1]+rate*(verts[i*9+7]-verts[i*9+1]); outfile<<mp[p]<<" "; p++; mp[p]=verts[i*9+2]+rate*(verts[i*9+8]-verts[i*9+2]); /*mp[p]=z;*/ outfile<<mp[p]<<"\n"; p++; num2+=2; continue; } } outfile.close(); batch.Begin(GL_LINES, num2); batch.CopyVertexData3f(mp); batch.End(); } void getstlmodel() { float boundary1[2]; float boundary2[2]; float boundary3[2]; bool isbegin=false; long size=0; int nlines=0; int count1=0; int count2=0; FILE* file=fopen("mystl.stl","r"); fseek(file,0L,SEEK_END); size=ftell(file); fclose(file); file=fopen("mystl.stl","r"); for (int i=0;i<size;i++) { if(getc(file)=='\n') { nlines++; } } num=nlines/7; rewind(file); while (getc(file) != '\n'); verts=new float[9*num]; vnorms=new float[9*num]; for (int i=0;i<num;i++) { char x[200]=""; char y[200]=""; char z[200]=""; if(3!=fscanf(file,"%*s %*s %80s %80s %80s\n",x,y,z)) { break; } vnorms[count1]=vnorms[count1+3]=vnorms[count1+6]=atof(x); count1++; vnorms[count1]=vnorms[count1+3]=vnorms[count1+6]=atof(y); count1++; vnorms[count1]=vnorms[count1+3]=vnorms[count1+6]=atof(z); count1+=7; fscanf(file,"%*s %*s"); if (3!=fscanf(file,"%*s %80s %80s %80s\n",x,y,z)) { break; } verts[count2]=atof(x); count2++; verts[count2]=atof(y); count2++; verts[count2]=atof(z); count2++; if (isbegin==false) { isbegin=true; boundary1[0]=boundary1[1]=atof(x); boundary2[0]=boundary2[1]=atof(y); boundary3[0]=boundary3[1]=atof(z); } else { if (verts[count2-3]>boundary1[0]) { boundary1[0]=verts[count2-3]; } else if (verts[count2-3]<boundary1[1]) { boundary1[1]=verts[count2-3]; } if (verts[count2-2]>boundary2[0]) { boundary2[0]=verts[count2-2]; } else if (verts[count2-2]<boundary2[1]) { boundary2[1]=verts[count2-2]; } if (verts[count2-1]>boundary3[0]) { boundary3[0]=verts[count2-1]; } else if (verts[count2-1]<boundary3[1]) { boundary3[1]=verts[count2-1]; } } if (3!=fscanf(file,"%*s %80s %80s %80s\n",x,y,z)) { break; } verts[count2]=atof(x); count2++; verts[count2]=atof(y); count2++; verts[count2]=atof(z); count2++; if (verts[count2-3]>boundary1[0]) { boundary1[0]=verts[count2-3]; } else if (verts[count2-3]<boundary1[1]) { boundary1[1]=verts[count2-3]; } if (verts[count2-2]>boundary2[0]) { boundary2[0]=verts[count2-2]; } else if (verts[count2-2]<boundary2[1]) { boundary2[1]=verts[count2-2]; } if (verts[count2-1]>boundary3[0]) { boundary3[0]=verts[count2-1]; } else if (verts[count2-1]<boundary3[1]) { boundary3[1]=verts[count2-1]; } if (3!=fscanf(file,"%*s %80s %80s %80s\n",x,y,z)) { break; } verts[count2]=atof(x); count2++; verts[count2]=atof(y); count2++; verts[count2]=atof(z); count2++; if (verts[count2-3]>boundary1[0]) { boundary1[0]=verts[count2-3]; } else if (verts[count2-3]<boundary1[1]) { boundary1[1]=verts[count2-3]; } if (verts[count2-2]>boundary2[0]) { boundary2[0]=verts[count2-2]; } else if (verts[count2-2]<boundary2[1]) { boundary2[1]=verts[count2-2]; } if (verts[count2-1]>boundary3[0]) { boundary3[0]=verts[count2-1]; } else if (verts[count2-1]<boundary3[1]) { boundary3[1]=verts[count2-1]; } fscanf(file,"%*s"); fscanf(file,"%*s"); } GLint xvert[108]= { 0,0,0, 1,1,0, 1,0,0, 1,1,0, 0,0,0, 0,1,0, 1,0,0, 1,0,1, 0,0,0, 1,0,1, 0,0,1, 0,0,0, 0,0,1, 1,0,1, 1,1,1, 0,0,1, 1,1,1, 0,1,1, 0,0,1, 0,1,1, 0,1,0, 0,0,1, 0,1,0, 0,0,0, 0,1,1, 1,1,1, 1,1,0, 0,1,1, 1,1,0, 0,1,0, 1,1,1, 1,0,1, 1,0,0, 1,1,1, 1,0,0, 1,1,0 }; int num1=0; for (int num2=0;num2<36;num2++) { xverts[num1]=boundary1[xvert[num2*3]]; num1++; xverts[num1]=boundary2[xvert[num2*3+1]]; num1++; xverts[num1]=boundary3[xvert[num2*3+2]]; num1++; } shaderManager.InitializeStockShaders(); boundaryBatch.Begin(GL_TRIANGLES, 36); boundaryBatch.CopyVertexData3f(xverts); boundaryBatch.CopyNormalDataf(xnormal); boundaryBatch.End(); std::cout<<"xmax="<<boundary1[0]<<"\txmin="<<boundary1[1]<<"\n"; std::cout<<"ymax="<<boundary2[0]<<"\tymin="<<boundary2[1]<<"\n"; std::cout<<"zmax="<<boundary3[0]<<"\tzmin="<<boundary3[1]<<"\n"; float step=(boundary3[0]-boundary3[1])/n; for (float i=boundary3[1]+step;i<boundary3[0];i+=step) { getpline(i,num,polyline[count],verts); count++; } } void ProcessMenu(int value) { switch(value) { case 1: tip=1; glutPostRedisplay(); break; case 2: tip=2; glutPostRedisplay(); break; case 3: tip=3; glutPostRedisplay(); break; case 4:; tip=4; break; } glutPostRedisplay(); } void SetupRC() { // Black background glClearColor(0.3f, 0.3f, 0.3f, 1.0f ); shaderManager.InitializeStockShaders(); viewFrame.MoveForward(1000.0f); triangleBatch.Begin(GL_TRIANGLES, num*3); triangleBatch.CopyVertexData3f(verts); triangleBatch.CopyNormalDataf(vnorms); triangleBatch.End(); boundaryBatch.Begin(GL_TRIANGLES, 36); boundaryBatch.CopyVertexData3f(xverts); for (int i=0;i<108;i++) { xverts[i]=-xverts[i]; } boundaryBatch.CopyNormalDataf(xnormal); boundaryBatch.End(); // Make the torus } void SpecialKeys(int key, int x, int y) { if(key == GLUT_KEY_UP) viewFrame.RotateWorld(m3dDegToRad(-5.0), 1.0f, 0.0f, 0.0f); if(key == GLUT_KEY_DOWN) viewFrame.RotateWorld(m3dDegToRad(5.0), 1.0f, 0.0f, 0.0f); if(key == GLUT_KEY_LEFT) viewFrame.RotateWorld(m3dDegToRad(-5.0), 0.0f, 1.0f, 0.0f); if(key == GLUT_KEY_RIGHT) viewFrame.RotateWorld(m3dDegToRad(5.0), 0.0f, 1.0f, 0.0f); // Refresh the Window glutPostRedisplay(); } void ChangeSize(int w, int h) { // Prevent a divide by zero if(h == 0) h = 1; // Set Viewport to window dimensions glViewport(0, 0, w, h); viewFrustum.SetPerspective(35.0f, float(w)/float(h), 1.0f, 2000.0f); projectionMatrix.LoadMatrix(viewFrustum.GetProjectionMatrix()); transformPipeline.SetMatrixStacks(modelViewMatix, projectionMatrix); } void RenderScene() { glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); glEnable(GL_CULL_FACE); glEnable(GL_DEPTH_TEST); modelViewMatix.PushMatrix(viewFrame); GLfloat vRed[] = { 1.0f, 0.0f, 0.0f, 1.0f }; //shaderManager.UseStockShader(GLT_SHADER_FLAT, transformPipeline.GetModelViewProjectionMatrix(), vRed); shaderManager.UseStockShader(GLT_SHADER_DEFAULT_LIGHT, transformPipeline.GetModelViewMatrix(), transformPipeline.GetProjectionMatrix(), vRed); switch(tip) { case 1: glPolygonMode(GL_FRONT_AND_BACK,GL_FILL); triangleBatch.Draw(); break; case 2: glPolygonMode(GL_FRONT_AND_BACK,GL_FILL); triangleBatch.Draw(); glPolygonMode(GL_FRONT_AND_BACK, GL_LINE); boundaryBatch.Draw(); break; case 3: glEnable(GL_LINE_SMOOTH); glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glPolygonMode(GL_FRONT_AND_BACK, GL_LINE); glLineWidth(2.5f); for (int i=0;i<count;i++) { polyline[i].Draw(); } break; } modelViewMatix.PopMatrix(); glutSwapBuffers(); } int main(int argc, char* argv[]) { std::cout<<"請輸入要剖分的份數\n"; std::cin>>n; gltSetWorkingDirectory(argv[0]); glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH | GLUT_STENCIL); glutInitWindowSize(800, 600); glutCreateWindow("Geometry Test Program"); glutReshapeFunc(ChangeSize); glutDisplayFunc(RenderScene); glutSpecialFunc(SpecialKeys); glutCreateMenu(ProcessMenu); glutAddMenuEntry("show stl model",1); glutAddMenuEntry("show boundary",2); glutAddMenuEntry("subdivision", 3); glutAddMenuEntry("Subdivision surface", 4); glutAttachMenu(GLUT_RIGHT_BUTTON); GLenum err = glewInit(); if (GLEW_OK != err) { fprintf(stderr, "GLEW Error: %s\n", glewGetErrorString(err)); return 1; } getstlmodel(); SetupRC(); glutMainLoop(); return 0; }
還望大家見諒。
示例圖片：
bearing實體模型
包圍盒
某條與Z法向平面交線
以下截取了某些數據點。僅僅是一小部分
39.9962 31.193 15.6757
40.431 31.1435 15.6757
39.9962 31.193 15.6757
39 31.3152 15.6757
40.431 31.1435 15.6757
41.3467 30.82 15.6757
42.5561 30.1707 15.6757
41.7464 30.6788 15.6757
41.7464 30.6788 15.6757
41.3467 30.82 15.6757
42.9772 29.8928 15.6757
43.1324 29.7641 15.6757
42.9772 29.8928 15.6757
42.9095 29.949 15.6757
43.1324 29.7641 15.6757
43.6553 29.2257 15.6757
43.6553 29.2257 15.6757
43.8835 28.9908 15.6757
44.0475 28.7774 15.6757
44.119 28.6843 15.6757
44.0475 28.7774 15.6757
43.8835 28.9908 15.6757。。。。。。。。。。
螺旋槳的實體模型
包圍盒
某條與Z法向平面交線