Code covered by the BSD License  

Highlights from
Mesh voxelisation

4.91667

4.9 | 14 ratings Rate this file 72 Downloads (last 30 days) File Size: 134 KB File ID: #27390
image thumbnail

Mesh voxelisation

by

 

28 Apr 2010 (Updated )

Voxelise a triangular-polygon mesh.

| Watch this File

File Information
Description

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)
 
 

Acknowledgements

This file inspired Intriangulation Which Points Are Inside A 3d Watertight Triangulation?.

MATLAB release MATLAB 7.10 (R2010a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (26)
21 Aug 2014 Ci  
28 Jul 2014 Pavel M.

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

14 Feb 2014 Christine

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

27 Nov 2013 John

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?

06 Sep 2013 Sven

@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 :)

08 Aug 2013 Jerome Crocco

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.

08 Aug 2013 Jerome Crocco

Programs like this make me love matlab even more.

Great work, there when I needed it.

Cheers,

06 Aug 2013 Abishek Balsamy Kamaraj  
21 Jun 2013 Anusha

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

08 Feb 2013 Yun-Ho Jang

Hello Adam,

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

08 Feb 2013 Adam A

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

07 Feb 2013 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.

02 Feb 2013 Yun-Ho Jang

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

07 Dec 2012 xianjiaoda  
16 Oct 2012 Adam A

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

09 Oct 2012 Eff

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

09 Jul 2012 Adam A

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.

05 Jul 2012 Ryan W

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?

27 Apr 2012 Adam A

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');

26 Apr 2012 Alexandre Zuquete Guarato

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!!!

31 Jan 2012 Prakash Manandhar  
02 Dec 2011 Alan AF  
23 Mar 2011 James

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!

23 Feb 2011 K B  
07 Feb 2011 Adam A

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);

07 Feb 2011 Amir Reza Riazi

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
28 Apr 2010

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

14 May 2010

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

20 May 2010

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

02 Jun 2010

Provided a clearer example of how to use the code.

14 Jun 2010

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

15 Jun 2010

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

22 Jun 2010

Improved the automatic detection of binary and ascii STL files.

23 Jun 2010

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

23 Jun 2010

Bug fix.

23 Jun 2010

And another minor bug fix.

28 Jun 2010

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

07 Sep 2010

The grid no longer has to entirely enclose the mesh.

26 Nov 2010

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

02 Dec 2010

Fixed bug in automatic grid generation.

07 Feb 2011

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)

03 Mar 2011

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

31 Aug 2011

Minor edits to the documentation

12 Feb 2013

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.

12 Feb 2013

Added a missing sub-function

Contact us