View License

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

» Watch video

Highlights from
Mesh voxelisation

4.9 | 23 ratings Rate this file 116 Downloads (last 30 days) File Size: 134 KB File ID: #27390 Version: 1.20
image thumbnail

Mesh voxelisation


Adam A (view profile)


28 Apr 2010 (Updated )

Voxelise a triangular-polygon mesh.

| Watch this File

File Information

Polygon mesh voxelisation

Adam H. Aitkenhead
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.

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

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.

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.

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

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

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


This file inspired Stl Tools and Intriangulation(Vertices,Faces,Testp,Heavytest).

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 (36)
30 Nov 2016 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.

04 Nov 2016 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.

01 Nov 2016 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

29 Jul 2016 good

good (view profile)

13 Jun 2016 Asif Arain

31 May 2016 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!

29 Apr 2016 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? :)


Comment only
31 May 2015 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!

28 May 2015 Reshma Bhat

11 Feb 2015 Sicong Wang

21 Aug 2014 Ci

Ci (view profile)

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

Comment only
14 Feb 2014 Christine

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

Comment only
27 Nov 2013 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?

Comment only
06 Sep 2013 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 :)

Comment only
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.


06 Aug 2013 Abishek Balsamy Kamaraj

21 Jun 2013 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

Comment only
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

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,

Comment only
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.

Comment only
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

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,

Comment only
09 Oct 2012 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:

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


Comment only
09 Jul 2012 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.

Comment only
05 Jul 2012 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?

27 Apr 2012 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');

Comment only
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

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!

23 Feb 2011 K B

K B (view profile)

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

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

Comment only
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.


28 Apr 2010 1.1

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

14 May 2010 1.2

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

20 May 2010 1.3

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

02 Jun 2010 1.4

Provided a clearer example of how to use the code.

14 Jun 2010 1.5

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

15 Jun 2010 1.6

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

22 Jun 2010 1.7

Improved the automatic detection of binary and ascii STL files.

23 Jun 2010 1.8

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

23 Jun 2010 1.9

Bug fix.

23 Jun 2010 1.10

And another minor bug fix.

28 Jun 2010 1.11

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

07 Sep 2010 1.12

The grid no longer has to entirely enclose the mesh.

26 Nov 2010 1.14

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

02 Dec 2010 1.15

Fixed bug in automatic grid generation.

07 Feb 2011 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)

03 Mar 2011 1.17

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

31 Aug 2011 1.18

Minor edits to the documentation

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

12 Feb 2013 1.20

Added a missing sub-function

Contact us