Surface Area at given hight of an STL file

Hello All,

I have an STL file and I would like to go through it (z-direction) calculating the surface area on the X-Y plane.

I would like to export the data of the layer number and the corrosponding area measured.

My knowlage in matlab is limited so any tips would help.

Thanks!

Answer 1

[DGM]

Here is what I put together. This should work for nonconvex objects, objects with holes, and it should handle cases where triangles are tangential to the slice plane. Tangent faces are included in the slice area (by default).

I haven't thoroughly tested it for accuracy, but it seems about right. Do your own testing. I've included some simple test models.

The given function stlslicer() takes a model or model data, optional parameters, and returns a polyshape() array. The model can be specified as numeric FV lists, a triangulation object, or as a full path to an STL file. Data is expected to be triangulated. Don't expect it to work with a quad mesh or anything with arbitrary length faces.

Section properties such as area, perimiter, etc. can be extracted from the polyshape array as needed. By default plotting is disabled. Obviously the animated plotting doesn't work on the forum, so only the last frame is shown.

% needed for the forum
unzip slicerdemo.zip

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% slice the model at one or more specific z-levels
PS = stlslicer('bwfrust.stl','zlevels',0,'doplot',true,'layout','horizontal');

% find the area of all the slices
A = area(PS)

A = 80.2827

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% slice a different model at specific z-levels
PS = stlslicer('saddleblock.stl','zlevels',[0 0.5 1 1.25],'doplot',true,'layout','horizontal');

% find the area of all the slices
A = area(PS)

A = 4×1

1.0000 0.9739 0.9622 0.2406

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% slice a different model, but just do uniformly-spaced slices
% note that the tangent faces are highlighted
% also note that stlslicer() can return zlevels (whether given or internally calculated)
[PS zlevels] = stlslicer('stepholecube.stl','numslices',5, ...
    'doplot',true,'layout','horizontal');

●

% find the area of all the slices
A = area(PS);
[zlevels A]

ans = 5×2

-0.5000 0.7500 -0.2500 0.7500 0 0.7500 0.2500 0.4375 0.5000 0.4375

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

If you don't need the display for anything, you can run without it. The plot routine otherwise wastes a lot of time and system resources.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% slice the same model without any display
[PS zlevels] = stlslicer('stepholecube.stl','numslices',5);

% find the area of all the slices
A = area(PS);
[zlevels A]

ans = 5×2

-0.5000 0.7500 -0.2500 0.7500 0 0.7500 0.2500 0.4375 0.5000 0.4375

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

If for some reason, you don't want to include tangent faces, you can specify that.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the same thing, but exclude tangent faces
[PS zlevels] = stlslicer('stepholecube.stl','numslices',5,'includetangent',false);

% find the area of all the slices
A = area(PS);
[zlevels A]

ans = 5×2

-0.5000 0 -0.2500 0.7500 0 0.4375 0.2500 0.4375 0.5000 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

EDIT:

Made some improvements to the tool after some showerthinking. Added the 'includetangent' option, fixed a few leftover bits of temporary patchwork I'd forgotten, and fixed a nasty unforeseen bug.

I decided to feed a pile of STLs to the slicer, testing at random slice heights, and I found that occasionally SurfaceIntersection() would spuriously return near-duplicate vertices only when queried at very particular zlevels (to within +-eps(zlevel)). The problem occurred at zlevels nowhere near model vertices, so it appears to simply be a rounding issue, not a problem resulting from vertex or edge alignments.

These duplicates would cause closed paths (which are index lists) to be left open and then subsequently discarded, resulting in wildly wrong area calculations. Finding near-duplicate vertices and repairing the candidate edge list costs some time, but it appears to be working consistently now.

These pitfalls just goes to show that this task wasn't exactly a simple problem that could be solved with some bullet points and a referral to some webdocs pages. Even relying on high-level tools like polyshape() for polygon simplification, this still amounts to over 1300 lines of code.

This uses this file from the FEX (included):

https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection

These are prior incomplete attempts to do the same:

https://www.mathworks.com/matlabcentral/answers/526274-how-to-get-the-cross-sectional-area-of-a-pde-model

https://www.mathworks.com/matlabcentral/answers/364768-how-to-post-process-output-data-of-surfaceintersection-function

[Juan]

This is pretty useful, thanks for sharing.

Answer 2

[Dev]

Hi @Omar,

We can follow the following steps to calculate the surface area in the X-Y plane, while traversing along the Z-direction-

• Read the STL file using the ‘stlread’ function in MATLAB. For more details on the usage of same, please refer to the link below- https://www.mathworks.com/help/releases/R2022b/matlab/ref/stlread.html
• The above function returns a triangulation object which can be used to extract the vertices and faces of the STL model read in the previous step. For more information on triangulation, please refer to the link below-https://www.mathworks.com/help/releases/R2022b/matlab/ref/triangulation.html
• Next, to traverse in the Z-direction (i.e., Z-coordinate), we can extract and store the unique Z-values to further loop over them. I have attached a sample code snippet below which achieves the same.

% Determine the unique Z-values
uniqueZ = unique(vertices(:, 3));

• For surface area calculation, we identify the faces that lie in the Z-plane since area calculation needs to be carried out in the Z-direction only. We can do this by looping through each unique Z-value and create a logical index for faces where all vertices are on the current Z-plane. We can then store these faces to finally use them for computing area.
• Next, to calculate the surface area, we can loop over the selected faces (identified in the previous step) and fetch the vertices of each selected face. Using these vertices, we compute area using the ‘polyarea’ function in MATLAB. For a detailed usage of this function, please refer the documentation below-https://www.mathworks.com/help/releases/R2022b/matlab/ref/polyarea.html
• Finally, we can store the layer number along with its corresponding area together in a 2D numeric array. I have attached a reference code snippet for the same below-

% Store the layer number and corresponding area
layerInfo(i, :) = [i, totalArea];

I hope the above steps help resolve the query.

[DGM]

The prior answer had convexity problems, but this answer doesn't even appear to be complete.

You propose to get a list of faces which have vertices belonging to a unique z-level defined strictly by the z-levels represented by the set of all vertices. That limits our z-sampling to however the object was triangulated. A simple object such as a quadrilateral prism constructed with 8 vertices and level end faces cannot be sampled at any intermediate z-level.

You then propose to use polyarea() to find the area represented by the vertices belonging to these selected faces. Polyarea() only accepts 2D inputs, so it's completely unclear how you're going to take a bunch of vertices (which aren't all at the same z-level) and find the area.

Consider the following simple saddle block (attached):

There are only three unique z-values in the set of all vertices. The only faces which have all three vertices at the same z-level are the two triangles on the bottom. How would you apply polyarea() to find the sectional area at z = 1? What about z = 1.25?