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<<"
";
    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]<<"
";
            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]<<"
";
            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]<<"
";
            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]<<"
";
            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]<<"
";
            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]<<"
";
            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)=='
')
        {
            nlines++;
        }
    }
    num=nlines/7;
    rewind(file);
    while (getc(file) != '
');
    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
",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
",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
",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
",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]<<"	xmin="<<boundary1[1]<<"
";
    std::cout<<"ymax="<<boundary2[0]<<"	ymin="<<boundary2[1]<<"
";
    std::cout<<"zmax="<<boundary3[0]<<"	zmin="<<boundary3[1]<<"
";
    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<<"请输入要剖分的份数
";
    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
", 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法向平面交线