arraysmatlab3dstl-format

Load stl file in matlab and convert to a 3D array


I have a stl file and I've loaded it in Matlab using stlread function. At this point I have a set of faces and vertices. How can I convert these faces and vertices in a 3D binary array like 512x512x100 array to obtain a binary 3D volume?


Solution

  • Ah lucky you. I am working with STL files recently and I coded some functions to do exactly this.

    First, note that you lose precision. STL files represent arbitrary shapes with arbitrary precision and converting it into a volume results in discretization and losses.

    That said, there is a very easy method to know if something is inside or outside a closed, connected triangulated surface, regardless if its convex or not: Throw a ray to the infinite and count intersection with the surface. If odd, its inside, if even, outside.

    The only special code you need is the line-triangle intersection, and the Möller Trumbore algorithm is one of the most common ones.

    function in=inmesh(fv,points)
    %INMESH tells you if a point is inside a closed,connected triangulated surface mesh
    % Author: Ander Biguri
    maxZ=max(fv.vertices(:,3));
    counts=zeros(size(points,1),1);
    for ii=1:size(points,1)
        ray=[points(ii,:);points(ii,1:2) maxZ+1];
        for jj=1:size(fv.faces,1)
            v=fv.vertices(fv.faces(jj,:),:);
            if all(v(:,3)<ray(1,3))
                continue;
            end
            isin=mollerTrumbore(ray, fv.vertices(fv.faces(jj,:),:));
            counts(ii)=counts(ii)+isin;
        end
    
    end
    in=mod(counts,2);
    end
    

    From FileExchange, with small modifications:

    function [flag, u, v, t] = mollerTrumbore (ray,tri)
    % Ray/triangle intersection using the algorithm proposed by Moller and Trumbore (1997).
    %
    % IMPORTANT NOTE: Assumes infinite legth rays.
    % Input:
    %    ray(1,:) : origin.
    %    d : direction.
    %    tri(1,:), tri(2,:), tri(3,:): vertices of the triangle.
    % Output:
    %    flag: (0) Reject, (1) Intersect.
    %    u,v: barycentric coordinates.
    %    t: distance from the ray origin.
    % Author: 
    %    Jesus Mena
    
        d=ray(2,:)-ray(1,:);
        epsilon = 0.00001;
    
        e1 = tri(2,:)-tri(1,:);
        e2 = tri(3,:)-tri(1,:);
        q  = cross(d,e2);
        a  = dot(e1,q); % determinant of the matrix M
    
        if (a>-epsilon && a<epsilon) 
            % the vector is parallel to the plane (the intersection is at infinity)
            [flag, u, v, t] = deal(0,0,0,0);
            return;
        end
    
        f = 1/a;
        s = ray(1,:)-tri(1,:);
        u = f*dot(s,q);
    
        if (u<0.0)
            % the intersection is outside of the triangle
            [flag, u, v, t] = deal(0,0,0,0);
            return;          
        end
    
        r = cross(s,e1);
        v = f*dot(d,r);
    
        if (v<0.0 || u+v>1.0)
            % the intersection is outside of the triangle
            [flag, u, v, t] = deal(0,0,0,0);
            return;
        end
        if nargout>3
            t = f*dot(e2,r); % verified! 
        end
        flag = 1;
        return
    end
    

    Just generate your points:

    yourboundaries=% get the range of your data from the STL file.
    [x,y,z]=meshgrid(yourboundaries);
    P=[x(:) y(:) z(:)];
    in=inmesh(fv,P);
    img=reshape(in,yourboundariesSize);