File Exchange

image thumbnail

Mesh voxelisation

version (134 KB) by Adam A
Voxelise a triangular-polygon mesh.


Updated 12 Feb 2013

View License

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)

Cite As

Adam A (2020). Mesh voxelisation (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (66)

Great contribution. But i have a question: how do i save the output to be used in another Simulator?
thanks in advance

Naveen B

Thank you sir it is useful for me Great contribution from your side

Hi, great app! it seems to work, but the results i get from it have many missing voxels. My 3d image has bubbles and i need to perform a morphological close operation to get a decent bubble. Also, i can't get W. Warriner's files to work. Any help is appreciated.

Yoeri Boink

I think I found another bug in the same function as in my previous comment. If an edge is vertical and a ray passes through that edge, then the side-check value is NaN because of a 0/0 division. The test is inadequate to catch these cases. It's possible these cases were caught by the interpolation error correction. That said, checking specifically for these cases allows for more robust single-ray or single-plane voxelization, allowing for faster and denser rasterized slicing of 3D geometries.

There may be a bug with respect to edge checking. If a ray arrives exactly at an edge, it is not considered part of the associated face. This can produce unexpected missing voxels in e.g. the "typical" triangulation of a cube centered at the origin with an odd number of equal-length voxels in each axial direction. In this case, one voxel is centered precisely at the origin. The rays sort of "sneak" between the triangles on every face, and that voxel gets missed. The error appears to be in lines 491, 497, and 503. The comparisons are strict (>). Changing them to (>=) appears to fix this issue. I am unsure if this will cause other unexpected issues, but it seems unlikely given how double precision arithmetic works. I fixed the issue in a github repo:



How can I save the resulting data to use as a Finite Element Mesh?

Miwa Sumiya


Thank you for this contribution. When I try to run : [gridOUTPUT] = VOXELISE(gridX,gridY,gridZ,STLin) on my mesh model having ~ 15 000 vertices and ~ 30 000 faces, I get a big Matlab memory error : Error using false
Requested 14961x14960x14960x3 (9355.0GB) array exceeds maximum array size preference. Creation
of arrays greater than this limit may take a long time and cause MATLAB to become
unresponsive. See array size limit or preference panel for more information.

Error in VOXELISE (line 281)
gridOUTPUT = false( voxcountX,voxcountY,voxcountZ,numel(raydirection) );

Error in compute_raw_skeleton (line 23)
vox_mesh =

Is this normal ? Is your code able to manage that number of vertices and faces ? (not that big to me). Or did I make something wrong ? Thank you for answer.


Huilin Yang

It's very helpful, thanks a lot.

Javed Was

Thank you for your work - I really appreciate it.

Hi Adam, thank you very much for this excellent tool. Pardon my possible ignorance, I am in need of extracting the coordinates of the full 3D voxelized mesh. In the trial of doing this, I presumed I should plot the gridCOx,y and z coordinates. Upon plotting this grid, instead of plotting the whole mesh, simply a line of points is plotted. Am I missing something?

hi adam,

this function has help me a lot. I just have a question, when I enter this command "[gridOUTPUT,gridCOx,gridCOy,gridCOz] = VOXELISE(gridX,gridY,gridZ,STLin,raydirection)" it gives me the grid x, y and z coordinate, but somehow it doesn't start at 0 or the origin.

Xu Chen

George Yang

Hi, it seems that after creating the voxel the scale has changed?

Hi Adam,
Thanks for the upload! Are there any papers/references on which the toolbox is based? Could you share them?



I also tried this and result I got,was completely wrong.
I tried with this one:
meshX = [3;0;0];
meshY = [0;3;0];
meshZ = [0;0;3];
[OUTPUTgrid] = VOXELISE(3,3,3,meshX,meshY,meshZ);

In the variable [Output grid],I could see only zeros :( :(.Plz help!


I tried to run the code using this:
meshX = [0;5;3];
meshY = [0;0;3];
meshZ = [0;0;0];
[OUTPUTgrid] = VOXELISE(6,6,6,meshX,meshY,meshZ);

I dont see any output here. I have a facet,which exists only in 2D Plane.Can you please help me to use your function


Maybe change line 301 in VOXELISE.m to:

gridOUTPUT = sum(gridOUTPUT,4)>=1;

This seems to have fixed my problem with corner-Voxels (which are lying on the corner of the grid) not beeing detected.

ntsh kr

how can i find the dimension of a voxel( i need to find out volume of a voxel)

Paul D



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


Asif Arain


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!


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


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


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


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


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?


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



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

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,

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!


Adam A

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


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?


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.

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?

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

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


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!


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.

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



Added a missing sub-function

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.

Minor edits to the documentation

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

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)

Fixed bug in automatic grid generation.

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

The grid no longer has to entirely enclose the mesh.

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

And another minor bug fix.

Bug fix.

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

Improved the automatic detection of binary and ascii STL files.

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

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

Provided a clearer example of how to use the code.

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

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

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

MATLAB Release Compatibility
Created with R2010a
Compatible with any release
Platform Compatibility
Windows macOS Linux