version 1.3.0.0 (6.39 KB) by
Peter Hammer

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

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.

Peter Hammer (2021). Marching Cubes (https://www.mathworks.com/matlabcentral/fileexchange/32506-marching-cubes), MATLAB Central File Exchange. Retrieved .

Created with
R2011a

Compatible with any release

- AI, Data Science, and Statistics > Statistics and Machine Learning Toolbox > Industrial Statistics >
- Industries > Medical Devices > DICOM Format >
- MATLAB > Graphics > 2-D and 3-D Plots > Surfaces, Volumes, and Polygons > Volume Visualization > Scalar Volume Data >
- Engineering > Biomedical Engineering > Biomedical Imaging >

**Inspired:**
BiofilmQ, resolution truncated with 10 digit'

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

张杰 张Leon CheungAram EskandariOk, I got it to work now. How does this algorithm solve ambiguities though, does it use maybe a midpoint/asymptotic decider?

Aram EskandariHey, I'm trying to use this function on a few regular cubes, but I can't get it to work. The first cube I have has values (see https://imgur.com/Mnv6UBs) [C, G, H, D] = [-1, 3, 2, 3] and [A, E, F, B] = [1, 2, -2, 2], and all my c-values are just 0. What I tried to do codewise is writing

C = [-1 3; 3 2];

C(:,:,2) = [1 2; 2 -2]; % creates a 3D matrix with two floors

[nY,nX,nZ] = size(C);

[X,Y,Z] = meshgrid(1:nX,1:nY,1:nZ);

[F,V,col] = MarchingCubes(X,Y,Z,C,0);

I don't get any errors from this, but the returning vectors F, V, col are all empty, and the resulting figure is also empty. Why is this, how am I using the algorithm wrong?

Peter HammerIn order for a reader to help, you must describe exactly how you called the function - including the dimensions of the inputs.

halilhi

I have the same problems as the one Muneeb that is

when i run this code then the folwing error is occur anybody help me to solve this error

Error using marching (line 33)

Manh DinhJIAWEI QINnaveen Bzelin1030hi Peter,

I don't know how to get the [x,y,z,c,iso] from point location data. Can you provide a method?

thanks

Morten Bormann NielsenA small bug I found: When you remove the duplicate vertices, you do not remove the same entries from the interpolated color list. This results in too many output values in the color array. A simple fix is to add "col = col(M);" immediately below "V = V(M, :);".

muneeb ullahhi

when i run this code then the folwing error is occur anybody help me to solve this error

Error using marching (line 33)

wrong number of input and/or output arguments

Hang Jo Annhi Peter,

I have 100 ultrasound images (cross section area of vessel).How can i apply your function for my project?Can you provide me some sample code? Thanks a lot!

Peter HammerAlex - the face list (F) and vertex list (V) output by marching cubes are, in general, not the same length, and they are not required to be the same length in order to construct a patch object from them using the patch function.

Pete

Alex MagsamPeter,

My output F and V are different lengths and cannot be used in the 'patch' function. Do you know why this would be?

In context, I have a 3D object I am rotating using roll, pitch and yaw transformations. The F and V are calculated correctly with the original coordinates, but after the transformation is applied F and V are outputted at different lengths.

Thanks,

Alex

yunbai wanghello,

THanks for the code and it worked great! But I have have a problem of the surface location.

The locations and the shape of surface dose not match with the original image volume.

I think the problem is that the voxel size of my image volume is not [1 1 1], it is [1.89 1.89 1].

how could I fixed this problem?

I really need to save the surface in the right position！

Many thanks

darshan ruikarhello,

really thank you very much for this code

i have a question

how to use COLORS input parameter and if we supplied that what i will get in cols

I have 586 CT images having communicated fractures . I implemented seed region growing algorithm on it to identify individual fractured zones out of it per fragment i assigned 1, 2, 3, . . . as fragment numbers

i want to color those fragment using different colors

how do i????

TikoTikoHello,

Thanks for the code. Can you please specify what data types are acceptable for matrix c of values? I am using this to visualize binarized data and have noticed there is a slight difference depending on whether my data is of type uint8 (0-background, 1-foreground) or of type logical, my ISO is 0.5. I have managed to trace this back to the line 181 of your code, where negative values will get truncated in case of uint8, where in case of logical they get converted to double. I imagine that the logical version is correct, as the intended data type according to your documentation is double or float.

Best regards!

Peter HammerHi Manika - You raise an important point. Matrices are not indexed in x,y,z order. The first and second indices of matrix C correspond to its rows and columns, respectively. To make coordinate vectors that are passed to meshgrid correspond to appropriate dimensions of matrix C, you must define the dimensions of matrix C this way:

[nY,nX,nZ] = size(C); % note first output arg corresponds to dimension in y-direction

Then simply call meshgrid and MarchingCubes exactly as you did below.

-Pete

ManikaManikaHi Peter,

my data set C has dimension 217*181*181

when I am writing:

1. [nX,nY,nZ] = size(C);

2. [X,Y,Z] = meshgrid(1:nX,1:nY,1:nZ);

3. [F,V,col] = MarchingCubes(X,Y,Z,C,4);

It is showing error saying X, Y, Z, and C should agree on dimensions...

but meshgrid by default generates X, Y Z with dimensions 181*217*181...

how to get C with 181*217*181 when original matrix dimensions are 217*181*181

Please help.

Thanks & Regards

zee falconI've download this project and try to run it in MATLAB 6.1 but it's not running(showing no output). Is there any thing else i've to download in matlab to run this??

Kindly guide me how can i run this code and see result??

Fatma KhansrcFiles = dir('E:\be project\BRAINIX\BRAINIX\IRM cérébrale, neuro-crâne\T1-3D-FFE-C - 801\*.dcm'); % the folder in which images exists

n = length(srcFiles);

a={}; slices={};

for i = 1:n

filename = strcat('E:\be project\BRAINIX\BRAINIX\IRM cérébrale, neuro-crâne\T1-3D-FFE-C - 801\',srcFiles(i).name);

a{i}= dicomread(filename);

slices{i} = (a{i});

end

X = cat(3,slices{1:end});

map = hsv(90);

XR = X;

Ds = smooth3(XR);

Sir, this is my code through which i have included the database in my project but i do not understand how to include marching cubes after that.

Alexander HeinzCorrection to the previous post:

Remove duplicate faces and keep order:

F = unique(F,'rows','stable');

Alexander HeinzF contains duplicate faces. Fixed thru:

% ---- Remove duplicate faces -------

FF = sortrows(F);

mask = [true; any(diff(FF),2)];

F = F(mask,:);

% -----------------------------------

Alexander HeinzThe faces contains duplicate vertices. The problem lies in the line 148:

F = I(F);

I fixed it provisionally by removing of faces with duplicates:

dupl_list = [];

for i=1:size(F,1)

uniq_faces = unique(F(i,:));

counter = histc(F(i,:),uniq_faces);

value = uniq_faces(counter > 1);

if ~isempty(value)

dupl_list = [dupl_list; i];

end

end

F = removerows(F,'ind',dupl_list);

Peter HammerHi Akash - I needed to do the same thing and found the following File Exchange contribution that worked great:

splitFV by Sven Holcombe (contributed in May 2010, updated June 2014)

Good luck!

-Pete

AKASH GHow to identify and separate the number of objects in the 3D mesh created?

MWasHi

You can remove the small holes in the surface by using the following function of removing duplicate vertices instead of the function (by Oliver Woodford) in marching cube code.

Just delete that section in code and call this function.

It works for me. Thanks

function [SV,SVI,SVJ] = remove_duplicate_vertices(V,epsilon,varargin)

% REMOVE_DUPLICATE_VERTICES Remove duplicate vertices upto a uniqueness

% tolerance (epsilon)

%

% SV = remove_duplicate_vertices(V,epsilon)

% [SV,SVI,SVJ] = ...

% remove_duplicate_vertices(V,epsilon,'ParameterName',ParameterValue,...)

%

% Inputs:

% V #V by dim list of vertex positions

% epsilon uniqueness tolerance (significant digit)

% Optional:

% 'WhiteList' Only merge vertices from the following selection (not

% working correctly, yet)

% Outputs:

% SV #SV by dim new list of vertex positions

% SVI #SV by 1 list of indices so SV = V(SVI,:)

% SVJ #V by 1 list of indices so V = SV(SVJ,:)

%

% Example:

% % Mesh in (V,F)

% [SV,SVI,SVJ] = remove_duplicate_vertices(V,1e-7);

% % remap faces

% SF = SVJ(F);

%

% default values

whitelist = [];

% Map of parameter names to variable names

params_to_variables = containers.Map( ...

{'WhiteList'}, ...

{'whitelist'});

v = 1;

while v <= numel(varargin)

param_name = varargin{v};

if isKey(params_to_variables,param_name)

assert(v+1<=numel(varargin));

v = v+1;

% Trick: use feval on anonymous function to use assignin to this workspace

feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));

else

error('Unsupported parameter: %s',varargin{v});

end

v=v+1;

end

if isempty(whitelist)

assert(nargin==1 || epsilon >= 0);

if nargin==1 || epsilon == 0

[SV,SVI,SVJ] = unique(V,'rows','stable');

else

[~,SVI,SVJ] = unique(round(V/(10*epsilon)),'rows','stable');

SV = V(SVI,:);

end

else

error('not implemented correctly')

VW = V(whitelist,:);

J = 1:size(V,1);

JW = J(whitelist);

[SVW,SVIW,SVJW] = remove_duplicate_vertices(VW,epsilon);

SJ = 1:size(V,1);

SJ(whitelist) = JW(SVJ);

end

end

BluppDear Peter, thank you for your work. However, I noticed that there still is an issue with duplicate vertices, which can create holes in the surface when vertices are moved.

Ninakavitha sunduHey, What i have is a three dimensional matrix with different grain id's . That is from 1 to 25 .The size of the matrix is 32*32*32. My question would how to adjust the iso value which satisfies all the grain id values. Is there any work around. Thanks a lot in advance. Looking forward to hear from you Peter.

Sidan YaoDear Peter, thank you so much for your work! It helped a lot!!

thirumagal jayaramanActually I have set of MRI images. 19 slices of the brain image. Using the marching cubes algorithm I created a 3D mesh by stacking all the slices. Now I want to get the histogram of this mesh having 19 slices-entire brain. Is it possible to get the histogram of the 19 slices using the mesh?

Gianni Schena..but if you include also the normals as additional output you deserve 6 stars

Athul sukumarHello Peter,

I want to convert 2d mri brain image to 3d image. how to use this code for executing the above problem?

Peter HammerDear Ahmed, I cannot understand your question. Is it related to my marching cubes function?

-Pete

Peter HammerDear Thirumagal, what exactly do you mean by plotting a histogram of the mesh? The mesh is a list of the Cartesian coordinates of points in 3D plus a list of triangles, identified by the indices of the points that comprise them. If a histogram is a set of sorted bins, what is the quantity that you want to sort?

-Pete

ahmed kareem1985Hi I need ,code to compute photo exterior orientation please

Parameshthirumagal jayaramanI want to plot histogram of the output mesh. But the output of Marching cubes is a list of faces and vertices. How can a store it in a variable and plot histogram?

Anna WhelanPeter HammerDear Sherly - I do not understand your question. Are you asking for the syntax of how to call the function with appropriate input and output arguments? Are you asking if I can reduce the entire algorithm to a "formula"? Please clarify.

-Pete

Sherly sariHello peter dou you have formula from your marching cubes function?

Peter HammerDear Jayant,

Marching cubes does not create a 3D mesh - rather it returns a 2D surface that intersects the 3D, volumetric data such that voxels on one side of the surface have an intensity value above a threshold value while those on the other side have an intensity below the threshold value. It does the same thing as Matlab's "isosurface" function. The primary advantage to my implementation of marching cubes over Matlab's isosurface is that the former is faster and the method is straightforward/transparent.

-Pete

jayant prajapatI have CT image(128 x 128 gray scale) of brain tumor. I want to create 3D mesh of the tumor. How to use this Marching Cube for this purpose

Sarah MillerHi Peter - Per your response to Eleo on 24 Jul 2015, if a binary stack (0s and 1s) is converted to a grayscale stack (0s and 255s), running the code with an iso of 127 (or anything but 0 or 255) should get you the correct results. Can you please confirm this? Thanks!

atmaxHi Peter,

Could you tell me how to calculate the number of cube containing surface and how to change the cube size please?

Best

Chris

RyanAwesome! Thanks so much!

-Ryan

Peter HammerRyan - If I understand your problem correctly, you could use the example code fragment below to create two separate surface models with marching cubes then plot them on the same axes in different colors.

Pete

[f1,v1] = MarchingCubes(X1,Y1,Z1,C1,ISO1);

[f2,v2] = MarchingCubes(X2,Y2,Z2,C2,ISO2);

figure

ph1=patch('vertices',v1,...

'faces',f1,...

'facecolor',[1 0 0],...

'facelighting','phong');

hold on

ph2=patch('vertices',v2,...

'faces',f2,...

'facecolor',[0 1 0],...

'facelighting','phong');

axis equal

light

RyanHi Peter,

I was interested to see if there is a way that I could plot multiple isosurfaces at once with different colors on the same plot?

Thanks,

-Ryan

Peter HammerAli - I am having a difficult time understanding exactly what you are trying to do. Marching cubes is a very straightforward algorithm: you pass it a 3D grayscale image along with the grayscale intensity value at which you want to construct a surface (i.e., the isovalue). The function returns the surface, defined as a list of triangles whose vertices are determined by interpolating between image voxels that lie on opposite sides of the isovalue. It does not select the isovalue automatically - you must specify it. Is this what you are trying to do?

Pete

ali HasanDear Peter

I am wondering how can I use marching cube algorithm for segmentation because I already use snake algorithm for segment tumors and I want to compare it with marching cube

Regards

Ali

Peter HammerDear Ali - I am glad you were able to use marching cubes with your data. To increase smoothness, you could simply insert a slice between every two adjacent slices of your original MRI data that is the mean of the two adjacent slices.

Pete

ali HasanDear Peter

I success to implement your code on segmented MRI slices and it gives me good 3D model but I need to make the 3D model more smooth because the margin of slices is appear clearly on 3D model do you know how can I make interpolation between slices in order to increase the smoothness between slices?

Regards

Ali

ali HasanDear Peter

Thank you for your reply

Is it possible to use Marching cube to segment the tumors such as other segmentation methods or not?

Regards

Ali

Peter HammerDear Ali,

No, marching cubes is not for building 3d models. It is for fitting a triangulated isosurface to 3d volumetric data. You can think of it as segmenting the 3d image...which it sounds like you have already done.

Pete

ali HasanDear Peter

Thank you for your code, I segment the tumor over a series of MRI slices successfully and I want to built a 3D model is it possible to use your code to build my 3D model or not

Regards

Ali

Peter HammerDear Eleo T,

Marching cubes operates on 3d grayscale images, not on segmented volumes. Try running it on the MRI data, with isovalue set to the same value used to arrive at your segmentation.

-Pete

Eléo TDear Peter,

Thank you for your function, it helps a lot !

I've to reconstruct three segmented volumes of a MRI one (WM, GM and CRL) (dimensions : 181 X 217 X 181)

I've used your function and tried a lot to work it out without success ...

I'm thinking I'm not far for it to work, but Matlab returns me this result :

F=[ ]

V=[ ]

I've put iso=1, (according to your last comment, I've done the his to and obtained 2 pics : one in 0 and one in 1, because my volumes are segmented from the original MRI volume).

Can you help ?

Thanks for your work

Eleo

Peter HammerDear Karen,

The MarchingCubes function is used simply by executing the single command:

[F,V] = MarchingCubes(x,y,z,c,iso);

Note that x, y, and z are matrices of the same size as c, which is the 3-d matrix of image intensity (grayscale) values. You can create coordinate matrices x, y, and z using Matlab's meshgrid command:

[x,y,z]=meshgrid(xmin:xdelta:xmax,ymin:ydelta:ymax,zmin:zdelta:zmax);

Here, xmin and xmax are the x-coordinates of the lower and upper extents of the 3-d image and xdelta is the x resolution. (Similarly for y and z). In the call to MarchingCubes, iso is the isovalue, i.e., the value of image intensity at which you want the output surface to lie. For example, for a CT scan in which you want to segment the surfaces of bone, iso is a value of image intensity that lies between the black background and the darkest elements that you consider to be bone. There is no universal method for selecting this value, as it depends on exactly what the user hopes to achieve with the segmentation. One suggestion is to look at the histogram of image intensity values with a command like this:

hist(c(:),200)

Sometimes the shape of the histogram shows 2 or more distinct peaks which are identifiable as, for example, bone and soft tissue. In this case just choose iso as a value that separates the two peaks.

KarenHello, I was wondering if you have an example script available somewhere that shows how to use this function correctly? (Along with sample DICOM files?) Thanks for sharing your work!

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

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

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

Peter HammerDear 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

Peter HammerDear 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

Marcos AlbarracinHi, 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

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

Peter HammerDear 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?

iulHello Dr,Peter,

Hello all,

Thank you very much for your code.

can u help me in my request of 3 jan 2015.

Regards;

Peter HammerDear 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

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

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

Peter HammerRao - 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?

nageswara raohi 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?

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

vivek KandasamyHi 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?

Itzik Ben ShabatPeter HammerVedpal - 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.

dai hongThanks Peter!

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

Peter HammerVedpal - 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

Vedpal SinghHello,

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

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

Itzik Ben ShabatIs 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.

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

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

VDThank you, Peter.

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

VDHi 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?

Anton SemechkoThanks for the submission Peter! It works great!

VDThank you very much Peter.

Peter HammerHi 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

VDHi,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?

Anton SemechkoI 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?

MarleenHi Peter,

Never mind the question below, I figured it out.

Thanks for your m-file, it works perfectly!

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

zhonglee323@gmail.comhello 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

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

Sarah MillerQuestion: 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.

ganghi,Peter,Tks for your work!

Is it the same as the funciton 'isosurface' of matlab 2012b?

Chen Victhanks a lot!

Oliver WoodfordBrilliant. Very fast.

Shusenwugood