How to make a STL file from a mesh

I am trying to make a STL file for a 3d printer. I have made 3d meshes of topographical data for various states and put nans in for points outside of the state boundary. I would like to convert these meshes (lonGrid, latGrid, El) into an STL. Additionally I would like to fill in all of the space below the mesh down to an elevation of 0.

I have tried http://www.mathworks.com/matlabcentral/fileexchange/42876-surf2solid-make-a-solid-volume-from-a-surface-for-3d-printing but it cannot seem to handle non-rectangular surfaces. My nans for out of state points break it.

I also tried http://www.mathworks.com/matlabcentral/fileexchange/4512-surf2stl but when I try loading it into Blender it does not look right. It just looks like random noise.

If I can get the mesh into a triangulated patch with the underside filled in I think I can get that converted to an STL.

Any help would be appreciated. Thanks

[Luis Jimenez]

Hey, a little late

Try this: stlwrite

Answer 1

[DGM]

Well. I thought for sure I'd used those before and it would be super easy to throw together a demo. I was wrong. As far as I can tell, this would have been quite a pain. I'm going to assume since these are geographic regions, there's no reason to presume that the boundary is convex. For now though, let's start with a simplified example. We want to create a round saddle shape.

% a surface
n = 30;
x = linspace(-1,1,n);
[X Y] = meshgrid(x);
Z = X.*Y + 1;

% truncate the surface with NaNs
r = hypot(X,Y);
Z(r>1) = NaN;

% observe the pringle
surf(X,Y,Z);
axis equal

The problem is as you say; sur2stl() doesn't extend the surface, and surf2solid() won't accept anything with NaNs in the vertex list. All we have is gridded data. We don't even have a triangulation, and triangulation of (potentially) non-convex regions is problematic. That said, we don't need to triangulate the region. Just do the entire mesh and delete the unwanted F,V data.

Once we have unambiguous FV data with no NaNs, surf2solid() works just fine. All that needs to be done is write it to an STL. In 2014, there would have been STL writing tools on the forum (e.g. Sven Holcombe's stlTools). Since R2018b, the base toolbox has stlwrite().

% triangulate the projection of the entire mesh
% this could be done with delanuay(), but it's overkill.
% FEX mesh2tri() was available, and it'll be faster.
% that's what surf2solid() would be using internally.
%V = [X(:) Y(:) Z(:)]; % if using delaunay()
%F = delaunay(V(:,1:2)); % this isn't really necessary
[F,V] = mesh2tri(X,Y,Z,'f'); % 'f','b','x';

% prune faces which have NaN vertices
% get rid of NaN vertices when done
goodverts = find(~isnan(V(:,3)));
goodfaces = all(ismember(F,goodverts),2);
F = F(goodfaces,:); % get rid of junk faces
map = zeros(size(V,1),1); % remap vertex indices
map(goodverts) = 1:numel(goodverts);
V = V(goodverts,:); % get rid of junk vertices
F = map(F); % remap faces to match new vertex list

% we can use surf2solid() now that it's unambiguously triangulated
[F,V] = surf2solid(F,V,'elevation',0);

% write it to a file
% stlWrite('pringle.stl',F,V) % from stlTools on FEX (before R2018b)
stlwrite(triangulation(F,V),'pringle.stl') % R2018b+

% display it using patch()
figure
patch('faces',F,'vertices',V, ...
    'facecolor','w','edgecolor','k');
view(3); camlight; view(10,33)
axis equal; grid on
xlabel('X'); ylabel('Y'); zlabel('Z')

Okay, so that's a start. So what's the problem? The problem is that stl2solid() presumes that the region has only one boundary. It will break and create a mess if there are holes or multiple unconnected regions. In this case, it creates a slot in the part and results in a bunch of invalid geometry.

It could be fixed internally, but at this point in the day, it's just easier to start over with a simplified example than it is to try figuring out someone else's code. Let's try a different simple example.

% a surface
n = 30;
x = linspace(-1,1,n);
[X Y] = meshgrid(x);
Z = X.*Y + 1;

% truncate the surface with NaNs
r = hypot(X,Y);
Z(r>1) = NaN;
Z(r<0.5) = NaN; % add a hole

% observe the holy pringle
surf(X,Y,Z);
axis equal

% triangulate the projection of the entire mesh
[F,V] = mesh2tri(X,Y,Z,'f'); % 'f','b','x'

% prune faces which have NaN vertices
% get rid of NaN vertices when done
goodverts = find(~isnan(V(:,3)));
goodfaces = all(ismember(F,goodverts),2);
F = F(goodfaces,:); % get rid of junk faces
map = zeros(size(V,1),1); % remap vertex indices
map(goodverts) = 1:numel(goodverts);
V = V(goodverts,:); % get rid of junk vertices
F = map(F); % remap faces to match new vertex list

% now the fun starts
% convert to a triangulation object so we can get 
% the open boundary as ordered lists
T = triangulation(F,V);
BV = freeBoundary(T); 

% test edges to see if they're unordered and reassemble them if necessary.
% if they're ordered, but concatenated, just split them.
sample = BV(1:end,:);
sample(:,2) = circshift(sample(:,2),1);
spans = diff([1; find(diff(sample(2:end,:),1,2)); size(BV,1)]);
isordered =  nnz(spans < 3)/numel(spans) < 0.5; % could be adjusted
if isordered
    boundaries = splitedges(BV); % attached
else
    boundaries = assembleedges(BV); % attached
end

% go through the boundaries and generate wall triangles
% we're relying on the fact that the boundary vertex order
% given by freeBoundary() is properly handed
nverts = size(V,1);
nfaces = size(F,1);
wallF = cell(numel(boundaries),1);
for k = 1:numel(boundaries)
    va = boundaries{k}(:,1); % top
    vb = va + nverts; % bottom
    vc = circshift(va,-1); % top
    vd = vc + nverts; % bottom
    wallF{k} = [vc vb vd; vb vc va]; % for forwardslash [t b b; b t t]
    %wallF{k} = [va vb vd; vd vc va]; % for backslash [t b b; b t t]
    % there's no need for a 'cross' option since all these quads are planar.
end
F = [F; cell2mat(wallF)]; % add them to the list

% create the bottom vertices
% this is a bit wasteful, since the interior points are coplanar
% but it's simpler than trying to retriangulate it
baseelevation = 0; % pick something
V = [V; V];
V(nverts+1:end,3) = baseelevation;

% create the faces for the bottom by just flipping the winding
F = [F; fliplr(F(1:nfaces,:))+nverts];

% write it to a file
% stlWrite('holypringle.stl',F,V) % from stlTools on FEX (before R2018b)
stlwrite(triangulation(F,V),'holypringle.stl') % R2018b+

% display it using patch()
figure
patch('faces',F,'vertices',V, ...
    'facecolor','w','edgecolor','k');
view(3); camlight; view(10,33)
axis equal; grid on
xlabel('X'); ylabel('Y'); zlabel('Z')

That looks like it'll work just fine, but there's one more problem that would have made this not work in 2014. At least in R2015b, the edges returned by freeBoundary() aren't sequentially sorted. Because of that, everything breaks.

The problem with it is that I don't actually know whether there are other conditions where this can happen (e.g. particular orientations of nonconvex boundaries), or when the behavior changed, so I'm going to need to do some more digging and come back to this.

EDIT: I couldn't find anything about this change in the old documentation, so I don't know when it happened. It's more verbose, but I added an automated test to see if the edge list is unsorted and handle it differently as needed. I guess it goes to show that "a bit more of a pain in 2014" was on point.

EDIT AGAIN: I got rid of delaunay() and cleaned up a few things. It's a bit faster, and a bit more readable. Using mesh2tri() is more than adequate, and it's faster and period-correct. I should probably put this effort into fixing surf2solid(), but I don't really feel like it.

[Christine Tobler]

Hi,

At a quick look internally, it seems freeBoundary starts to have ordered outputs in R2018b.

[DGM]

Oh. I didn't actually dig into the files. I was digging through my stash of old RN.

Well. It's good to have that nailed down. Thank you.