File Exchange

image thumbnail

Mesh voxelisation

version 1.20 (134 KB) by

Voxelise a triangular-polygon mesh.

81 Downloads

Updated

View License

Polygon mesh voxelisation
=========================

Adam H. Aitkenhead
adam.aitkenhead@christie.nhs.uk
The Christie NHS Foundation Trust

Voxelize a closed (ie. watertight) triangular-polygon mesh. The mesh can be in one of several formats: in an STL file; in a structure containing the faces and vertices data; in three 3xN arrays containing the x,y,z coordinates; or in a single Nx3x3 array defining the vertex coordinates for each of the N facets.
 
 
USAGE:
======

[gridOUTPUT,gridCOx,gridCOy,gridCOz] = VOXELISE(gridX,gridY,gridZ,STLin,raydirection)
  ...or...
[gridOUTPUT,gridCOx,gridCOy,gridCOz] = VOXELISE(gridX,gridY,gridZ,meshFV,raydirection)
  ...or...
[gridOUTPUT,gridCOx,gridCOy,gridCOz] = VOXELISE(gridX,gridY,gridZ,meshX,meshY,meshZ,raydirection)
  ...or...
[gridOUTPUT,gridCOx,gridCOy,gridCOz] = VOXELISE(gridX,gridY,gridZ,meshXYZ,raydirection)
 
 
INPUT PARAMETERS
================

gridX - Mandatory - 1xP array - List of the grid X coordinates.
                    OR an integer - Number of voxels in the grid in the X direction.

gridY - Mandatory - 1xQ array - List of the grid Y coordinates.
                    OR an integer - Number of voxels in the grid in the Y direction.

gridZ - Mandatory - 1xR array - List of the grid Z coordinates.
                    OR an integer - Number of voxels in the grid in the Z direction.

STLin - Optional - string - Filename of the STL file.

meshFV - Optional - structure - Structure containing the faces and vertices of the mesh, in the same format as that produced by the isosurface command.

meshX - Optional - 3xN array - List of the mesh X coordinates for the 3 vertices of each of the N facets
meshY - Optional - 3xN array - List of the mesh Y coordinates for the 3 vertices of each of the N facets
meshZ - Optional - 3xN array - List of the mesh Z coordinates for the 3 vertices of each of the N facets

meshXYZ - Optional - Nx3x3 array - The vertex positions for each facet, with:
                    1 row for each facet
                    3 columns for the x,y,z coordinates
                    3 pages for the three vertices

raydirection - Optional - String - Defines in which directions the ray-tracing is performed. The default is 'xyz', which traces in the x,y,z directions and combines the results.
 
 
OUTPUT PARAMETERS
=================

gridOUTPUT - Mandatory - PxQxR logical array - Voxelised data (1=>Inside the mesh, 0=>Outside the mesh)

gridCOx - Optional - 1xP array - List of the grid X coordinates.
gridCOy - Optional - 1xQ array - List of the grid Y coordinates.
gridCOz - Optional - 1xR array - List of the grid Z coordinates.
 
 
EXAMPLES
========

To voxelise an STL file:
>> [gridOUTPUT] = VOXELISE(gridX,gridY,gridZ,STLin)

To voxelise a mesh defined by a structure containing the faces and vertices:
>> [gridOUTPUT] = VOXELISE(gridX,gridY,gridZ,meshFV)

To voxelise a mesh where the x,y,z coordinates are defined by three 3xN arrays:
>> [gridOUTPUT] = VOXELISE(gridX,gridY,gridZ,meshX,meshY,meshZ)

To voxelise a mesh defined by a single Nx3x3 array:
>> [gridOUTPUT] = VOXELISE(gridX,gridY,gridZ,meshXYZ)

To also output the lists of X,Y,Z coordinates:
>> [gridOUTPUT,gridCOx,gridCOy,gridCOz] = VOXELISE(gridX,gridY,gridZ,STLin)

To use ray-tracing in only the z-direction:
>> [gridOUTPUT] = VOXELISE(gridX,gridY,gridZ,STLin,'z')
 
 
NOTES
=====

- Defining raydirection='xyz' means that the mesh is ray-traced in each of the x,y,z directions, with the overall result being a combination of the result from each direction. This gives the most reliable result at the expense of computation time.
- Tracing in only one direction (eg. raydirection='z') is faster, but can potentially lead to artefacts if rays exactly cross a facet edge.
  
 
REFERENCES
==========

- This code uses a ray intersection method similar to that described by:
  Patil S and Ravi B. Voxel-based representation, display and thickness analysis of intricate shapes. Ninth International Conference on Computer Aided Design and Computer Graphics (CAD/CG 2005)
 
 

Comments and Ratings (40)

paulygon

Jae-Wook Nam

Adrian Becker

Rik Wisselink

Thank you very much for this code. It really helped me a lot.

I came across a weird thing, which I want to let you know about.
I have delineations in .wrl format that are read in to a triangulated mesh. The z-coordinates of the vertices are all discrete multiples apart, as they follow the slice thickness from the CT. When I then use this script, it may happen that slices are returned empty when the values in gridZ happen to be equal to one of the z-values. This is not always the case.
The fix I used is simply using a tiny offset in the z-direction, but the artifact took me by surprise.

Michael C.

I've used this function to voxelise an stl.file.

How would I go about using the outputs for a slice function in matlab? Is this even possible?

Please Help.

Shen DU

Hi,I'd like to say it's really an excellent work. And it helps me a lot to solve my case.
I want to know after I get the OUTPUTgrid, how can I access to the normal vector of the voxel in the surface?

Thank you so so so much

good

good (view profile)

Asif Arain

Andre

Andre (view profile)

First of all, thanks for this submission, as it proved very useful for my cause.

The function worked perfectly. The only issue I have is that the voxelised matrix that is being outputted in gridOUTPUT is flipped along the X-dimension. Can anyone suggest why this might be happening please?

The code I used follows:
% Voxelise STLmodels (faces and vertices)
fv = stlread(data.stlList(l).name);
    
% Obtain point cloud from vertex data
data.stlPtCld{l} = pointCloud(fv.vertices);
    
% Get dimensions of STL models
sz(1) = ceil(data.stlPtCld{l}.XLimits(2) - data.stlPtCld{l}.XLimits(1));
sz(2) = ceil(data.stlPtCld{l}.YLimits(2) - data.stlPtCld{l}.YLimits(1));
sz(3) = ceil(data.stlPtCld{l}.ZLimits(2) - data.stlPtCld{l}.ZLimits(1));
    
% Voxelise STL model
data.stlMask{l} = VOXELISE(sz(1), sz(2), sz(3), fv, 'xyz');
    
% Flip Voxelised stlMask since output from |VOXELISE| is flipped along the X-dimension
data.stlMask{l} = flip(data.stlMask{l}, 1);
    
Many thanks in advance for any help!

Andre

Andre (view profile)

I am getting the following error when trying to voxelise:

Operands to the || and && operators must be convertible to logical scalar values.

Error in VOXELISE (line 239)
if ~isempty(strfind(raydirection,'x')) && (min(gridCOx)>meshXmin || max(gridCOx)<meshXmax)

any help please? :)

Thanks!

Andrew Hundt

This one worked very well for reading the STL file. Other STL reading matlab scripts were either limited by binary or ascii format, or simply didn't work. Thanks!

Reshma Bhat

Sicong Wang

Ci

Ci (view profile)

Pavel M.

How to display the resultant [OUTPUTgrid] = VOXELISE(grid_s,grid_s,grid_s,'sample.stl','xyz'); in 3D?(Volumetric plot)...

Christine

Has anyone attempted using this function to somehow return the volume of the voxelized solid?

John

John (view profile)

MATLAB newbie here...

I have a STL file, and after importing it into ML, it dumps data into a structure with faces and vertices.

What I'm having problems with is how to use your function with the faces and vertices data I have.

Any tips?

Sven

Sven (view profile)

@Jerome, if you'll indulge a little cross-promotion... check out inpolyhedron()... it performs a similar job with a different approach, and I've been looking for some real-world examples with a heavy memory footprint to test it against. It might (no promises!) do well for your input.

@Adam, nice work :)

Nice work. happy to use this.

Took about an hour to do a 3200 x 400 x 400 array. some memory problems at times.

Overall works really well.

Programs like this make me love matlab even more.

Great work, there when I needed it.

Cheers,

Anusha

Anusha (view profile)

Nay I know how can I create the voxelisation using 3D points(x,y,z) of on the surface only. The 3D points should be the vertices. But how can I get the normals. Hope to get some help

Yun-Ho Jang

Hello Adam,

Thanks a lot. I can get a large size binary mask based on your comments.

Adam A

Adam A (view profile)

Hi Yun-Ho Jang,

Yes, the code can be used to get a binary mask at a specified z-depth by specifying the required z-depth as the input gridZ. This works as long as gridZ is not an integer value. (By integer I do not mean in terms of datatype (int16, int32, int64, etc), I just mean in terms of being a whole number.)

However, if gridZ is an integer value then the code presumes that it defines the number of voxels in the grid in the z direction. To work around this, you could voxelise at two z-depths defined in gridZ (as a two-element vector), then discard the slice you don't want to leave you with a single binary layer. Not an elegant way to do it unfortunately, but this wasn't an anticipated use-case for the code when it was designed.

Hope this helps,
Adam

Yun-Ho Jang

Hello Adam,

Your code is great, but when I tried to voxelize with a bit large number of grids, (eg. 1000x1000x100), it generated out of memory error. Is there any way to get a binary mask at specific z-depth?

Thank you in advance.

Yun-Ho Jang

I was looking for a way to convert 3d stl data into a volume and this function is great!

xianjiaoda

Adam A

Adam A (view profile)

Hi Eff,
This code is designed for triangular-polygon meshes only, but your example mesh consists of quadrilateral faces.
Hope this helps,
Adam

Eff

Eff (view profile)

Hi Adam,

I'd like to use your code for the creation of a volume from an isosurface defined by the vertices and the faces as used by patch. However, the VOXELISE method does not return the expected results. I have created a simple example in order to reproduce the problem:
http://pastebin.com/PkWRuyyZ

Am I doing something wrong or is this a bug in the code?

Thanks,
Eff

Adam A

Adam A (view profile)

I suspect the best way will be to rotate the STL before voxelisation, like you suggest. If you try to do it post-voxelisation by rotating the 3D grid then you'll need to do a lot of interpolating, which will be slow and could introduce errors or artifacts. As for parallel processing - I don't see any reason why it shouldn't work, but I haven't looked into it.

Ryan W

Ryan W (view profile)

It is a great program to convert 3D stl to voxels.
Is there a efficient way to rotate the voxel model for arbitrary angle? One way I think of is that I rotate the 3D stl model (calculate the new coordinates of each vertices of the facets) and voxelise the rotated 3D model. But the 'voxelise' itself is quite slow. Is there any way to rotate the voxel model itself? Will parallel processing possible in voxelise?

Adam A

Adam A (view profile)

Hi Alexandre. If you want to obtain only the surface voxels rather than the entire volume, try the following line after doing the voxelisation. Hope this helps. Adam
gridSHELL = gridOUTPUT & imdilate(1-gridOUTPUT,ones(3,3,3),'same');

I really liked your mesh voxelisation code. I would like to know how can I insert into gridOUTPUT all surface voxels instead of just inside mesh voxels. thanks!!!

Alan AF

James

James (view profile)

Great voxelising script, I've found it very useful. One thing I'd like to see if possible - if I try to voxelise multiple bodies that overlap, where they overlap the script returns 'void' voxels rather than 'solid'. Preferably I'd like to see these returned as solid as well. Cheers!

K B

K B (view profile)

Adam A

Adam A (view profile)

Hi Amir,
I've just uploaded a new version of the code, which will give a list of unique vertex coordinates using the following two lines of code. Hopefully it will be available for download shortly.
Adam

[coordVERTICES] = READ_stl(‘filename.stl');
[faces,vertices] = CONVERT_meshformat(coordVERTICES);

Hi Adam, how could I just have the xyz coordinates of the points of the triangles without duplication?

I just want the points for easier analysis.

thanks

Updates

1.20

Added a missing sub-function

1.19

Added checking of ray/vertex intersections, which reduces artefacts in situations where the mesh vertices are located directly on ray paths in the voxelisation grid.

1.18

Minor edits to the documentation

1.17

Improved method of finding which mesh facets can possibly be crossed by each ray. Up to 80% reduction in run-time.

1.16

Include the utility CONVERT_meshformat.m, which enables the vertex coordinates to be converted from an Nx3x3 array to faces,vertices info. Ie.
[faces,vertices] = CONVERT_meshformat(coordVERTICES)

1.15

Fixed bug in automatic grid generation.

1.14

Changed handling of automatic grid generation to reduce chance of artefacts.

1.12

The grid no longer has to entirely enclose the mesh.

1.11

Allow input to be a structure containing [faces,vertices] data, similar to the type of structure output by isosurface.

1.10

And another minor bug fix.

1.9

Bug fix.

1.8

Enable ray-tracing in any combination of the x,y,z directions.

1.7

Improved the automatic detection of binary and ascii STL files.

1.6

Fixed a bug introduced by the change to a logical output array.

1.5

Output is now a logical array, improving the meory efficiency.

1.4

Provided a clearer example of how to use the code.

1.3

Now optionally output the grid x,y,z coordinates. Robustness also improved.

1.2

Now also works with non-STL input. Changes also provide a significant speed improvement.

1.1

Fixed a bug which caused artifacts to appear for some STLs.

MATLAB Release
MATLAB 7.10 (R2010a)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video