STL模型讀取的續——STL模型與面的交線

xds1224@163.com

前些天，我一個交大研究生同窗給了我一個大做業，讓我嘗試下，內容大體以下：

可以讀取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法向平面交線