.stl 3D模型文件的读取计算，方法和程序实现（matlab版&&C++版&&python版）

0. 背景描述

3D模型、3D打印中很常见的一种文件格式.STL文件，其描述的主要就是其表面点所组成的三角面片的点坐标信息（vertex），和法向量（normal）。 如果单纯查看的话，很多软件都可以打开，诸如SolidWorks、3DMax，win10自带3d查看器、MITK等。win10自带的3D查看器，较为推荐。而我最早一般是用基于vtk的医学可视化软件，进行读取查看。

如果想要访问获取.STL文件中的点坐标数据，来进行后续更为自由的计算和处理，最好是需要程序实现的。其中须知，.STL文件一般有两种文件模式，即ASCII文本型和Binary二进制形式。其中ASCII形式的文件，很多文本查看软件便可打开，如记事本、notepad++、sublime等等。推荐notepad++，其大文件的打开速度很优秀。

想对这两种文件模式下的STL文件进行统一的读取计算，读取stl,返回表面点坐标（可附加计算内部点坐标）。程序操作如下：

1. C++版

基于vtk进行C++的简单实现即可，vtk的大牛们已经写好了vtkSTLReader的类，程序实现demo见另外的link:URL.

2. Python版

同样类似于C++版，还是基于vtk实现，程序实现demo见另外的link:URL.

3. Matlab版

MATLAB的版本程序最麻烦，详见下面的stlGeneralAccess函数，支持二进制和文本型两类STL的访问和计算，并且帮项目中的同门实现了计算内部点坐标的程序部分（方法参考了程序注释中的REFS.2），如果没有这个需要的，把对应的代码(计算inNodes部分)删掉或者注释掉即可。

NOTE: 程序逻辑也比较简单，就是读入STL文件，判断是ASCII还是Binary，然后分别处理，先获得三角面片的关系和点坐标，然后针对点坐标进行去重整理，得到纯粹的表面点即边界点的坐标，之后利用文献2（REFS.2）的方法计算得到内部点。

function [x, y, z, bdNodes, inNodes] = stlGeneralAccess(filename, dimSize)
% This function reads an STL file in binary format or ascii format.
% It returns the coordinates of both boundary nodes and interior nodes
% of the 3D polyhedron.
%
% partial REFs:
% 1. stlread() by Doron Harlev for binary STL file via mathswork.com
% 2. https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html
%
% Examples:
%        close all
%        clear all
%        clc

%        [x, y, z, bdNodes, inNodes] = stlGeneralAccess('./Region1.stl');
%        if 0
%           patch(x, y, z);
%        else
%           scatter3(bdNodes(:,1), bdNodes(:,2), bdNodes(:,3), 'MarkerFaceColor',[.75 .75 .0]);
%           hold on
%           scatter3(inNodes(:,1), inNodes(:,2), inNodes(:,3),'MarkerFaceColor',[.5 .0 .0]);
%        end
%
% Where
%       filename  ---   Full path of the input STL file.
%       dimSize   ---   [3*1]. The dimSize of coordinate, covering the 3D polyhedron.
%       x         ---   [3*numFacet]. X Position of three points in a triangular facet.
%       y         ---   [3*numFacet]. Y Position of three points in a triangular facet.
%       z         ---   [3*numFacet]. Z Position of three points in a triangular facet.
%       bdNodes   ---   [num_bdNodes*3]. The coordinates of boundary nodes.
%       inNodes   ---   [num_inNodes*3]. The coordinates of interior nodes.
%
% FUNCTION stlGeneralAccess. Version 6.1 Written by JeffZhang. AUG,2020.
% E-mail: jfsufeng@gmail.com || jfzhang2018@zju.edu.cn
% Dsp.CN:该函数支持STL文件的二进制形式和文本形式的通用访问和计算，可计算得到边界点和内部点

if nargin < 1
    error('Too few input arguments')
end
if nargout > 5
    error('Too many output arguments')
end
fid=fopen(filename);
if fid == -1
    error('File could not be opened, pls check the path.')
end
disp('>>>stlGeneralAccess.[START]')
%tic;toc;
nline = 0;
binary = -1;
ascii = -1;
binaryOrText = -1;
stlLines={};
while ~feof(fid)
    tline = fgetl(fid);
    nline = nline + 1;
    if contains(tline, 'endfacet')
        ascii = 1;
        binary = 0;
    end
    if nline >= 50
        break;
    end
end

if ascii == 1 && binary ==0
    disp('FileMode: text');
    binaryOrText = 0;
else
    disp('FileMode: binary');
    binaryOrText = 1;
end

num_facet = 0;
frewind(fid)
switch binaryOrText
    case 0
        while ~feof(fid)
            stlLines=[stlLines;fgets(fid)];
        end
        
        %norm = [];
        ver1=[];ver2=[];ver3=[];
        x=[];y=[];z=[];
        nmRow=[];
        %x=zeros(3,num_facet); y=zeros(3,num_facet); z=zeros(3,num_facet);
        for i = 1:size(stlLines,1)
            %tline = fgetl(fid);
            tl = stlLines{i};
            id = strfind(tl,'facet normal');
            if ~isempty(id)
            %if contains(tl,'facet normal')
                nmRow=[nmRow;i];
                num_facet = num_facet+1;
                %tk = tl(id+13:end); %15 or id+13(better&more stable)
                %nm = strsplit(tk,' ');
                %norm = [norm; str2double(nm)]; %% skip
                
                if i+2 <= size(stlLines,1) && i+3 <= size(stlLines,1) && i+4 <= size(stlLines,1)
                    ver1Str = stlLines{i+2};
                    ver2Str = stlLines{i+3};
                    ver3Str = stlLines{i+4};
                    id1 = strfind(ver1Str,'vertex');
                    id2 = strfind(ver2Str,'vertex');
                    id3 = strfind(ver3Str,'vertex');
                    if isempty(id1) || isempty(id2) || isempty(id3)
                    %if ~contains(ver1Str,'vertex')||~contains(ver2Str,'vertex')||~contains(ver3Str,'vertex')
                        error('Exception0x0815 of vertex.')
                    end                              
                    tk1 = ver1Str(id1+7:end);  %11 or id1+7(better&more stable)
                    tk2 = ver2Str(id2+7:end);
                    tk3 = ver3Str(id3+7:end);
                    nm1 = strsplit(tk1,' ');
                    nm2 = strsplit(tk2,' ');
                    nm3 = strsplit(tk3,' ');
                    ver1 = [ver1; str2double(nm1)];
                    ver2 = [ver2; str2double(nm2)];
                    ver3 = [ver3; str2double(nm3)];
                    x=[x;ver1(end,1) ver2(end,1) ver3(end,1)];  %% For patch validation
                    y=[y;ver1(end,2) ver2(end,2) ver3(end,2)];
                    z=[z;ver1(end,3) ver2(end,3) ver3(end,3)];
                end
            end
        end
        x=x';y=y';z=z';
        
    case 1
        
        ftitle=fread(fid,80,'uchar=>schar'); % Read file title
        num_facet=fread(fid,1,'int32'); % Read number of Facets
        
        fprintf('\nTitle: %s\n', char(ftitle'));
        fprintf('Num Facets: %d\n', num_facet);
        x=zeros(3,num_facet); y=zeros(3,num_facet); z=zeros(3,num_facet);
        ver1=zeros(num_facet,3); ver2=zeros(num_facet,3); ver3=zeros(num_facet,3);
        for i=1:num_facet
            norm=fread(fid,3,'float32'); % normal coordinates, ignored for now
            ver1(i,:)=fread(fid,3,'float32'); % vertex 1
            ver2(i,:)=fread(fid,3,'float32'); % vertex 2
            ver3(i,:)=fread(fid,3,'float32'); % vertex 3
            col=fread(fid,1,'uint16'); % color bytes
            
            x(:,i)=[ver1(i,1); ver2(i,1); ver3(i,1)]; % convert to matlab "patch" compatible format
            y(:,i)=[ver1(i,2); ver2(i,2); ver3(i,2)];
            z(:,i)=[ver1(i,3); ver2(i,3); ver3(i,3)];
        end
        
end
fclose(fid);

bdVertexes=[ver1;ver2;ver3];
bdNodes=unique(bdVertexes, 'rows');

%% 缺少图像Size，Spacing等信息,所以这里只是做个临时demo,届时传入实际的Size[3]和Spacing[3]即可如果posCrd坐标计算的话
% 如果直接只用idxCrd坐标的话,直接完美使用,不用传参.
imSize = [512 512 300];
imSpacing = [1 1 1]; % [.5 0.5 0.5]
origin = [0. 0. 0.];
if nargin < 2
    disp('Default dimSize below.')
    dimSize = imSize
end

minX_bd = min(bdNodes(:,1));
maxX_bd = max(bdNodes(:,1));
minY_bd = min(bdNodes(:,2));
maxY_bd = max(bdNodes(:,2));
minZ_bd = min(bdNodes(:,3));
maxZ_bd = max(bdNodes(:,3));
% posCrd = origin + idxCrd*imSpacing, for medical image computing;
% 如果用posCrd来计算;则步长h(i)=imSpacing(i)；rangeSz(i) = imSize(i)*imSpacing(i);
inNodesPre=[];
inNodes=[];
for i= 1:1:dimSize(1)
    for j= 1:1:dimSize(2)
        for k= 1:1:dimSize(3)
            tmpX = origin(1) + i * imSpacing(1);
            tmpY = origin(2) + j * imSpacing(2);
            tmpZ = origin(3) + k * imSpacing(3);
            if tmpX >= minX_bd && tmpX <= maxX_bd && ...
                    tmpY >= minY_bd && tmpY <= maxY_bd && ...
                    tmpZ >= minZ_bd && tmpZ <= maxZ_bd
                inNodesPre = [inNodesPre; tmpX tmpY tmpZ];
            end
        end
    end
end

if isempty(inNodesPre)
    fprintf('\n!ERROR!Warning:The dimSize may be not suitable.\n')
end

%% fix x，evaluate y,z. (pnpoly method)
for i= 1:1:size(inNodesPre,1)
    polyEdge2dId = find((bdNodes(:,1)>= inNodesPre(i, 1)-.5)&(bdNodes(:,1)<= inNodesPre(i, 1)+.4));
    vertY = bdNodes(polyEdge2dId,2);
    vertZ = bdNodes(polyEdge2dId,3);
    testY = inNodesPre(i,2);
    testZ = inNodesPre(i,3);
    c = 0;
    
    for k=1:1:length(vertY)
        if k == 1
            j = length(vertY);
        else
            j = k - 1;
        end
        if ( ((vertZ(k)>testZ) ~= (vertZ(j)>testZ)) && ...
                (testY < (vertY(j)-vertY(k)) * (testZ-vertZ(k)) / (vertZ(j)-vertZ(k)) + vertY(k)) )
            c = ~c;
        end
    end
    if ~mod(c,2)
        %disp('even -- outside')
    else
        %disp('uneven -- inside')
        inNodes = [inNodes; inNodesPre(i,:)];
    end
end
disp('>>>stlGeneralAccess.[OVER]')

end

最后附两个该函数调用的demo效果：

 

 

 

此文暂结~ 后续的更新可以到我的mathworks-FileExchange去check&update.

URL1: https://www.mathworks.com/matlabcentral/fileexchange/79191-stlgeneralaccess

URL2: https://github.com/JeffJFZ/stlGeneralAccess