4.83333
4.8 | 14 ratings Rate this file 98 Downloads (last 30 days) File Size: 6.39 KB File ID: #32506 Version: 1.3
image thumbnail

Marching Cubes

by

Peter Hammer (view profile)

 

11 Aug 2011 (Updated )

Use vectorized marching cubes algorithm to compute triangulated mesh of an isosurface from 3D matrix

| Watch this File

File Information
Description

This function uses a vectorized version of the marching cubes algorithm to compute a triangulated mesh of the isosurface within a given 3D matrix of scalar values at a given isosurface value. The output is a triangulated mesh specified in terms of a face list and a vertex list. The orientation of the triangles is chosen such that the normals point from the higher values to the lower values. Optional arguments COLORS ans COLS can be used to produce interpolated mesh face colors.

This function was used to create a surface mesh from the CT scan dataset for the Stanford bunny, a 461 x 339 x 330 array of floats. See uploaded image. Elapsed time on an AMD Opteron 64-bit computer with 4 GB of RAM is 24.7 seconds. For comparison, a surface mesh was computed from the same dataset using Matlab's isosurface function, and the elapsed time was 98.6 seconds.

This function was adapted for Matlab by Peter Hammer in 2011 by making minor syntax changes to an Octave function written by Martin Helm <martin@mhelm.de> in 2009 (http://www.mhelm.de/octave/m/marching_cube.m)

Revised 30 September, 2011 to add code by Oliver Woodford for removing
duplicate vertices.

MATLAB release MATLAB 7.12 (R2011a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (44)
30 Mar 2015 Peter Hammer

Peter Hammer (view profile)

Dear iul,
In my experience, the dicom format for CT images is not always the same. Matlab has some tools for reading dicom images; if they don't work, you might have to use matlab's fread. Good luck.

Comment only
29 Mar 2015 iul

iul (view profile)

Hi,Dr.Peter,
i succesfully use your code for the CT head data set from this site(http://graphics.stanford.edu/data/voldata/),meanwhile my goal is to use the code in a dicom file format image each image have a size of 515kb,how to load these data in order to made the 3d matrix.

Thanks for your great support.

11 Mar 2015 brhm

brhm (view profile)

I fixed the problem Mr. Hammer. Ty for answering. Btw I'm dealing with organ segmentation, that's why I had to convert my CT data to the form of a logical 3D matrix in order to use that logical matrix as the input of another algorithm. Also I have some questions about landmarking of a volume (liver, kidney, etc.). If you have a background, I'd like to send a detailed e-mail that includes a couple of images attached.

Thank you very much again.

Regards..

Comment only
06 Mar 2015 Peter Hammer

Peter Hammer (view profile)

Dear Marcos,
I am not sure what you mean by a "three-dimensional mesh". Marching cubes performs the same task as Matlab's isosurface function: it creates a surface mesh of triangles that passes through volumetric grayscale data at a constant grayscale value. If this is what you want, marching cubes might be an appropriate tool. If you are looking for a three-dimensional volumetric mesh (e.g., a connected mesh of tetrahedra that occupies volume), it is not an appropriate tool. Before working on your 2nd problem (example of passing the arguments to the code), it makes sense to ensure that marching cubes is appropriate for your task.
Pete

Comment only
06 Mar 2015 Peter Hammer

Peter Hammer (view profile)

Dear brhm, the marching cubes method was developed for segmenting volumetric grayscale data. For grayscale data, ISO defined the grayscale value at which to construct the isosurface. For logical data, verying the isovalue between zero and one will not have much effect on the ouput surface...because the isosurface will always pass through the same voxel edges. Why is your CT data in the form of a logical 3D matrix?
Pete

Comment only
06 Mar 2015 Marcos Albarracin

Hi, Dr. Peter,

Thank you very much for your code.
I have two problem.
Firts i am new in the topic.
My goal is to build a three-dimensional mesh from a set of tomographic images of a porous medium.

I would like to know how your code would help me in my goal.

On the other hand like you could give a small example of how to run your code because I'm not getting.

Very thanks

Comment only
25 Feb 2015 brhm

brhm (view profile)

Hello Dr. Peter,

Thank you very much for your code. I'm planning to use your code to get faces and vertices informations which are critical for me. I executed your code with different ISO values (between 0 and 1 ) by holding other parameters the same for my logical 3D matrix which is loaded by CT slices and couldn't figure out how it works. What's the deal with ISO? I'd be really grateful if you could explain me how ISO affects the output. Thank you for your helps in advance.

Regards..

[F,V] = MarchingCubes(X,Y,Z,C,???ISO???)

Comment only
12 Jan 2015 Peter Hammer

Peter Hammer (view profile)

Dear iul,
You must first load the CT slices into a single 3d matrix, much like Hwathanie has done in the code fragment in the comments below. Have you been able to do this?

Comment only
10 Jan 2015 iul

iul (view profile)

Hello Dr,Peter,
Hello all,
Thank you very much for your code.
can u help me in my request of 3 jan 2015.

Regards;

Comment only
05 Jan 2015 Peter Hammer

Peter Hammer (view profile)

Dear Hwathanie,
For the CT data that I read, it was necessary to use big-endian byte ordering when reading the file. E.G.,

data=fread(fid,[512,512],'uint16','ieee-be');

Does this solve your problem?

Sincerely,
Pete

Comment only
03 Jan 2015 iul

iul (view profile)

Hello Dr.Petter,

Could you help me plz.
i have multiple Ct-scan slices of paratoid galnde(dicom file),and i want to make a 3D reconstruction .so how to create my 3d grid from these slices???

Regards;

03 Jan 2015 Hwathanie

Hello Peter, Thank you very much for your code. I was trying to use it with the bunny data set provided in [http://graphics.stanford.edu/data/voldata/] with the following Matlab code.

for i = 1:361
fid = fopen(strcat('bunny\',num2str(i)));
v(:,:,i) = reshape(fread(fid, inf, '*uint16'), [512 512]);
fclose(fid);
end

[x y z] = meshgrid(1:512, 1:512, 1:361);
x = single(x);
y = single(y);
z = single(z);

[F,V] = MarchingCubes(x,y,z,v,500);

This took a long time and gave me an output like this [http://tinypic.com/r/22axb7/8]. Could you please tell how to correctly use your code.

09 Jul 2014 Peter Hammer

Peter Hammer (view profile)

Rao - The function is not flexible with input and output arguments. To properly choose the input arguments, please read the long section of commented text at the beginning of the mfile. I cannot answer your question about antenna design without knowing more about your approach. What type of image/data are you passing to marching cubes?

Comment only
08 Jul 2014 nageswara rao

hi peter,
it shows an error
Error using ==> MarchingCubes at 33
wrong number of input and/or output arguments
can i use it for antenna designs?

Comment only
07 Jul 2014 Peter Hammer

Peter Hammer (view profile)

Vivek - It could be that the isovalue needs to be adjusted. I have noticed that when image is somewhat noisy and I choose too low an isovalue, the resulting surface is noisy as you describe.

Comment only
06 Jul 2014 vivek Kandasamy

Hi Peter, I have tried to create a mesh from volumetric CT data. But the final surface has spiky triangles on the surface creating an uneven surface. How can i avoid this?

22 Jun 2014 Itzik Ben Shabat  
12 Mar 2014 Peter Hammer

Peter Hammer (view profile)

Vedpal - Marching cubes is an algorithm to produce a mesh of a surface at a given isovalue contained within 3d volumetric data. In other words, it identifies only bounding surfaces and does not fill the volume within the boundaries with "3-d" mesh elements like tetrahedra or cubes. You will have to come up with a meshing algorithm to do this...or find a software package that can do this.

Comment only
01 Mar 2014 dai hong

Thanks Peter!

27 Feb 2014 Vedpal Singh

Thanks Peter, for your valuable suggestion.

Regarding the 2nd question:
Using marching cube, I got a 3D result. But this result is not perfect because resultant 3D area is not perfectly filled. Boundary is clear but some inside space is vacant.
So, how we fill out this vacant space inside the 3D image boundary.

26 Feb 2014 Peter Hammer

Peter Hammer (view profile)

Vedpal - The construction of the edgeTable and triTable values is tricky. I think the original paper on marching cubes by Lorensen & Cline does a pretty good job of explaining. I will email you the paper. Regarding your second question on filling holes, I am not sure I understand what you are asking. Can you be more specific?
-Pete

Comment only
22 Feb 2014 Vedpal Singh

Hello,
Really a very good work.
Actually,I am confused in function

[edgeTable, triTable] = GetTables();

How can we calculate the edgeTable and triTable values.

and

How we can feel the holes in 3D desined through marching cube algorithm.

Thanks

31 Oct 2013 Peter Hammer

Peter Hammer (view profile)

Itzik - Yes, the OpenGL renderer runs faster than default zbuffer - possibly a lot faster depending on hardware. I use a Phong lighting model in MarchingCubes, and Matlab's OpenGL renderer does not support Phong. Using the OpenGL renderer with an alternate lighting model, like Gouraud, is a good choice for large meshes where speed is an issue - and effects of the lighting model are less noticeable.

Comment only
31 Oct 2013 Itzik Ben Shabat

Is there a particular reason why you didnt use the opengl renderer property for the figure? it runs much faster for me when adding it instead of the default renderer.

Comment only
21 Oct 2013 Peter Hammer

Peter Hammer (view profile)

JK - Marching cubes works on volumetric data, not on point clouds. If your point cloud data represents points on a surface, you might explore convex hull based triagulation methods. See Matlab's DELAUNAY function, for example.

Comment only
21 Oct 2013 JK Hwang

Hi Peter, I have tried to create mesh from 3D point cloud. How can I use your code to get my goal?

13 Aug 2013 VD

VD (view profile)

Thank you, Peter.

Comment only
12 Aug 2013 Peter Hammer

Peter Hammer (view profile)

VD - Rather than converting your (grayscale) CT data to binary, pass the grayscale data to the function. Marching cubes interpolates between nodal values of image intensity to produce a smooth isosurface.

Comment only
10 Aug 2013 VD

VD (view profile)

Hi Peter, I have database of about 100 CT scan images and I've converted it into binary images having bones as white. Using your package I am generating 3D skull from binary slices, however it shows step-wise(spiky) surface. looks like each slice as separate step. How can I get smooth continuous surface using this package?

Comment only
25 Jul 2013 Anton Semechko

Thanks for the submission Peter! It works great!

23 Jul 2013 VD

VD (view profile)

Thank you very much Peter.

Comment only
23 Jul 2013 Peter Hammer

Hi VD - Input arguments x, y, and z are expected to be the same as would be passed to Matlab's isosurface function. In fact, I posted the vectorized marching cubes function to provide a fast alternative to isosurface. The help section for the isosurface function gives a pretty clear explanation of what is expected for these input arguments: the arrays X, Y, and Z represent a Cartesian, axis-aligned grid. C contains the corresponding values at these grid points. The coordinate arrays (X, Y, and Z) must be monotonic and conform to the format produced by meshgrid. C must be a 3D volume array of the same size as X, Y, and Z.

I hope this is helpful.

Pete

Comment only
22 Jul 2013 VD

VD (view profile)

Hi,I have difficulty understanding what I need to pass as x,y,z & c matrices? I have 99 slices of binary image for which I need to generate surface. Would you please give me some example?

Comment only
14 Jun 2013 Anton Semechko

I think it may be helpful, for some people who are concerned about the efficiency of the code, if you added to your description the actual size of the CT dataset that was used to test the function. In other words, how big was the CT image of the bunny?

Also did you have a chance to benchmark this code against the built-in 'isosurface' function?

Comment only
10 Jun 2013 Marleen

Hi Peter,

Never mind the question below, I figured it out.
Thanks for your m-file, it works perfectly!

10 Jun 2013 Marleen

Hi Peter,

Thanks for the marching cubes algorithm!
I am not so familiar with imaging yet, and I wonder how I should create the matrix C? I know for each point in the mesh (x,y,z) whether it is inside the structure or not. How do I translate this to C, or where can I find more information on such a matrix?

Thanks a lot!

Comment only
09 Apr 2013 zhonglee323@gmail.com

hello Hammer,

what is "3D matrix of scalar values C" actually mean? If I get a sample point matrix, which ndims = 2, how can I transform them to 3D matrix?

Thanks

Comment only
09 Apr 2013 Peter Hammer

Peter Hammer (view profile)

Sarah - the algorithm assumes no knowledge of the data outside of the 3D matrix, so it does not "close" the surface on the boundary. For example, if you take a CT scan of a basketball and apply marching cubes to a subvolume containing only a hemisphere, it will return two disconnected surfaces: one hemispherical shell for the outside surface of the ball and one hemispherical shell for the inside surface of the ball.

Comment only
08 Apr 2013 Sarah Miller

Question: If I have an item that does not terminate at the edges of my 3D matrix (the scan is only a part of the sample), will this algorithm count the exterior of the matrix as well? Let me know if more info is needed.

Comment only
26 Jan 2013 gang

gang (view profile)

hi,Peter,Tks for your work!
Is it the same as the funciton 'isosurface' of matlab 2012b?

14 Nov 2012 Chen Vic

thanks a lot!

Comment only
30 Sep 2011 Oliver Woodford

Brilliant. Very fast.

01 Sep 2011 Shusen

Shusen (view profile)

 
20 Aug 2011 wu

wu (view profile)

good

Comment only
Updates
18 Aug 2011 1.1

Added image of Stanford bunny meshed using this function. Added text describing this example to description.

30 Sep 2011 1.2

Revised 30 September, 2011 to add code by Oliver Woodford for removing duplicate vertices.

23 Jul 2013 1.3

Changed description to include benchmarking results, and changed help section of mfile to more explicitly describe input arguments.

Contact us