Delaunay Triangulation and a Plane

When reconstructing a surface especially in three dimensional, its almost inevitable that you need to generate planes based on a set of vertices. In the case, when the facets are unknown, we basically have a cloud of three dimensional points. Based on the cloud of points, our eyes may be able to tell us what the sillouette or the contour of the three dimensional model in question is; however it is quite difficult or this is not trivial and it is a well-studied problem in computational geometry.

Certainly, there are a lot of recently proposed and well-developed methodologies for estimating or reconstructing a surface based on a given cloud of points. This post is written just to provide an idea of one of the basics: delaunay triangulation.  In short delaunay triangulation is to create a triangular mesh for a set of vertices that meets the Delaunay conditions(the main idea is to prevent any surface-within-a-surface).

In matlab, it is fair easy to follow the example (provided by Matlab).  Or if you want to explore further, you can try to make a cube. Basically you just need to write down the vertices and let Matlab do its job. I will leave the readers to do this on your own.

X, Y, Z are the column vectors of a cube(meaning all the x coords, y coords and z coords, say you have 2 vertices (0,0,0) and (0,0,1), X=[0 0]; Y=[0 0]; Z=[0 1];).

> dt = delaunay(X,Y,Z); 
> trimesh(dt,X,Y,Z); 

Now, that's pretty straightforward. However, if you have a problem similar to mine, when you only want to create a 3D plane with delaunay(you have your own reasons, i don't have to tell you mine). In the above example, if you only specify one of the faces of the cube, Matlab will probably tell you that they cannot proceed because all vertices are on a plane or are co-linear.

In order to do so, you need to trick Matlab into thinking the plane is no longer a plane; basically you just have to 'extrapolate' the plane in the perpendicular direction

Originally 
X =[0 1 1 0];
Y =[0 1 1 0];
Z =[0 0 1 1];

Perturbed 
X =[ 0 1 1.01 1.01 1 0];
Y =[0 1 0.99 0.99 1 0];
Z =[0 0 0 1 1 1];

Running the same command in the previous section should give you the following rendering.

Rendering a 3D Volume from a stack of 2D images

I have seen a number of questions related to the title.  I even experienced the question myself.  Here, i would like to document my approach.

The Question


If you are in the medical imaging field, at the time of writing, you are almost guaranteed that you'd need to reconstruct a 3D volume based on a stack of 2D images(let it be a stack of 'gray-scale' images).  The question is how in Matlab?


The Possible Solutions

---

Part 0 - Assumptions

---

• all 2D images are of the same dimension, hence your 3D volume can hold all of them in a rectangular cube
• majority of the pixels in each of the 2D images have 3D spatial relationships (you can't visualize much if the pixels in each of the 2D images are of some random distribution. )

---

Part 1 - Visualizing 3D Volume from A Stack of 2D Images

---

To visualize or reconstruct a 3D volume from a stack of 2D images, you can try the following toolkits in matlab.

[1] 3D CT/MRI images interactive sliding viewerhttp://www.mathworks.com/matlabcentral/fileexchange/29134-3d-ctmri-images-interactive-sliding-viewer

[2] Viewer3D http://www.mathworks.com/matlabcentral/fileexchange/21993-viewer3d

[3] Image3 http://www.mathworks.com/matlabcentral/fileexchange/21881-image3

[4] Surface2Volume http://www.mathworks.com/matlabcentral/fileexchange/8772-surface2volume

[5] SliceOMatic http://www.mathworks.com/matlabcentral/fileexchange/764

Note that if you are familiar with VTK, you can try this: [6] matVTKhttp://www.cir.meduniwien.ac.at/matvtk/

I am currently sticking with [5] SliceOMatic for its simplicity and ease of use. However, by default, rendering 3D is quite slow in Matlab. Turning on openGL would give faster rendering. (http://www.mathworks.com/help/techdoc/ref/opengl.html) Or simply put, set(gcf, 'Renderer', 'OpenGL').

---

Part 2 - Interpolating pixels in between the slices

---

To interpolate pixels in between the slices, you need to specify an interpolation method (some of the above toolkits have this capability / flexibility. Otherwise, to give you a head start, some examples for interpolation are bicubic, spline, polynomial etc..(you can work this out by looking up on google or google/scholar for interpolation methods much more specific to your problem domain).

---

Part 3 - 3D Pre-processing

---

Looking at the procedures, one may process the volumetric data by processing each of the 2D images first. In many advanced algorithms, or in true 3D processing, what you can do is to process the volumetric data in 3D domain first (simply put, you take the 26 neighbors or more in to account first.). Once this step is done, you can simply output the volumetric data into a stack of 2D images for cross-sectional viewing or supply to one of the aforementioned toolkits for 3D viewing or output to third party 3D viewing applications.

---

I have followed the above concepts for my own medical imaging research projects.

The above writing is documented here but was aiming at answering this question on StackOverflow.

Running Matlab codes on multicores

Running Matlab codes with multicores has never been this easy!

For your great, great, great work, my sincere thank you to you, Markus Buehren. Thank you.

[Multicore - Parallel Processing on Multiple Cores]

http://www.mathworks.com/matlabcentral/fileexchange/13775

To run it, it's fairly simple.

Say you are looping a function called DoSomething(x, y, z) 20 times.

In your standard procedures, you would have:
for i =1:20,

%% derivation/computation of x,y,z

DoSomething(x,y,z);

end

Now you will need to rewrite your codes, for example:

paramCell = cell(20,3);

for i = 1:20

paramCell{1, i} = x;%% of course you could merge the 3 args together!!

paramCell{2, i} = y;%%

paramCell{3, i} = z;%%

end

startmulticoremaster(@DoSomething, paramCell, settings); %% for settings just copy settings from the demo!

good, nifty script!!!! Thanks once again!