Hi all,

I want to get a skeleton of a 3D volume. I'm pretty sure I can do this in a round-about way, but I think there might be a shortcut.

The 3D volume I have is an STL file (think points connected in triangles or, say, the output from the isosurface() function) of a bifurcating tree (ie, imagine a tree trunk with some connected branches).

Here's one way that I think I can do:
1. Place a volume of points over my STL mesh.
2. Convert to a voxel volume (InPolyhedron on the fex will do this)
3. Run a 3D skeletonization on the voxel volume (http://www.mathworks.com/matlabcentral/fileexchange/4917-skeleton-in-3d will do this)

HOWEVER, I have this nagging suspicion that I can do it more simply *without* converting to voxels, and simply using the triangular mesh directly. My suspicion is that I can do some fancy work utilising MATLAB's DelaunayTri class along with the voronoiDiagram method. Reading through the help I can't quite tell if I'm on the right track or chasing down a red herring. Can someone more familiar with TriReps tell me:
A) Is it possible to do what I want?
B) What steps should I take?

For those that like test data to start from:
load tetmesh
trep = TriRep(tet, X);
[tri Xb] = freeBoundary(trep);
myFV = struct('vertices',xf,'faces',tri); % This is my STL input
figure, patch(myFV,'EdgeColor','k','FaceColor','r','FaceAlpha',.5), axis image, view(3), camlight

Any ideas about how to reduce a volume like this to its basic skeleton? (ie, a representation with 4 lines running down each bar, and a (near to a) circle at the top and the base?

Thanks,
Sven.

Sven:
Thanks for the test data. But it bombed:
??? Undefined function or variable 'xf'.
Error in =3D=3D> test at 16

Anyway, here is some code to get the medial surface of a 3D volume.
The algorithm is due to Steve Eddins of the imaging group at the
Mathworks:
I start with a 3D image called image3D:

		% Start the timer.
		[rows cols slices] =3D size(image3D);
		tStart =3D tic;
		% Compute the eroded 3D image because sometimes the particles touch.
		SE =3D strel('ball', 1, 1, 0);
		eroded3DImage =3D imerode(image3D, SE);

		% Convert the volume so that the air=3D1 and the particles=3D0.
% 		image3D =3D 1 - image3D;

		% Compute the 3D skeleton using ITK's routine.
		% Commented out - didn't seem to work as expected.
% 		disp('Starting generation of 3D skeleton...');
% 		image3DSkeleton =3D ITKSkeleton(image3D);

		%------------------------
		% Compute the 3D skiz using Steve Eddins proposed method.
		% There is a way of using a 3-D watershed transform that I think
matches closely your inflating balls metaphor.
		% Specifically:
		% 1.       Compute the distance transform of your 3-D blobs image
using bwdist.
		%          (Fortunately, bwdist runs much faster and takes up less
memory for the multidimensional case in R2009b than in earlier
releases.)
		% 2.       Compute the watershed transform of the distance transform
using watershed.
		%          The zero-valued pixels in the output of watershed will be
at the interfaces between the inflated balls.
		% I think this is what morphology people call the =93SKIZ,=94 or
=93skeleton by influence zones.=94  The interior of each inflated ball is
an =93influence zone.=94
		% Okay, so here we go......
		message =3D sprintf('Generating 3D distance transform on %d rows by %d
columns by %d slices 3D image...', rows, cols, slices);
		disp(message);
		set(handles.txtInfo, 'string', message);
		try
			%---- Compute the 3D Euclidean Distance Transform.
			try
				memory
				edtImage3D =3D bwdist(eroded3DImage); % Calculate EDT
				clear('eroded3DImage'); % Clear eroded3DImage to free up memory.
			catch ME
				errorMessage =3D sprintf('Error running bwdist()\nError Message:\n
%s', ME.message);
				disp(errorMessage);
				uiwait(warndlg(errorMessage, 'bwdist error', 'error'));
			end
			elapsedSeconds =3D toc(tStart);

			%---- Compute the 3D watershed.
			message =3D sprintf('\n%s took %.1f seconds\nGenerating 3D watershed
transform on %d rows by %d columns by %d slices 3D image...', message,
elapsedSeconds, rows, cols, slices);
			disp(message);
			set(handles.txtInfo, 'string', message);
			try
				memory
				labeledImage =3D watershed_exp(edtImage3D, 6); % Calculate the
watershed.
% 				clear('edtImage3D'); % Clear edtImage3D to free up memory.
			catch ME
				errorMessage =3D sprintf('Error running watershed()\nError Message:\n
%s', ME.message);
				disp(errorMessage);
				uiwait(warndlg(errorMessage, 'watershed error', 'error'));
			end
			elapsedSeconds =3D toc(tStart);

			%---- Get the pixels with a value of zero
			message =3D sprintf('\n%s %.1f total seconds elapsed so far\nFinding
zero-value voxels on %d rows by %d columns by %d slices 3D image...',
message, elapsedSeconds, rows, cols, slices);
			disp(message);
			image3DSkeleton =3D (labeledImage =3D=3D 0);
			clear('labeledImage'); % Clear labeledImage to free up memory.


As you can see, I originally tried a version from ITK but that didn't
work.  And take note that I called it the "medial surface" and not the
"skeleton."  Because in 3D the "skeleton" is not necessarily a thin
single pixel wide line - it can be a surface.  Imagine the space
between two flat walls.  If that were thinned down, you'd have a sheet
(surface) lying between the two walls, not a single line.

Steve's algorithm works well in my application where I wanted the
medial surface of the void space going between particles that are
suspended in a 3D volume.
Good luck,
-ImageAnalyst

Thanks ImageAnalyst.

Sorry about the bug. It was a typo that somehow always comes around when you paste in MATLAB code to a browser and then make the tiniest of changes :)  See below which is fixed.

My question was actually about a triangulated surface, but I am interested in your implementation for a 3D image input. As far as I can tell, the program you gave can be summarised as follows:
1. Calculate the bwdist() of your binary 3D image.
2. Create a labelled 3D watershed image
3. The ridges of the watershed (where label==0) is your skeleton.

I understand these steps in general, but I found that for the life of me, calling watershed on my distance image simply returns a label matrix full of ones. In your code, you're calling "watershed_exp". Is this a file exchange file, or should the standard IPT "watershed" function work just fine?

Thanks,
Sven.


% CREATE TRIANGULATED INPUT
load tetmesh
trep = TriRep(tet, X);
[tri Xb] = freeBoundary(trep);
myFV = struct('vertices',xf,'faces',tri); % This is my STL input
figure, patch(myFV,'EdgeColor','k','FaceColor','r','FaceAlpha',.5), axis image, view(3), camlight 

% CONVERT TO 3D BINARY IMAGE
bboxPts = [min(myFV.vertices,[],1);max(myFV.vertices,[],1)];
bboxDiff = ceil(diff(bboxPts));
xVec = linspace(bboxPts(1,1),bboxPts(2,1), 2*bboxDiff(1));
yVec = linspace(bboxPts(1,2),bboxPts(2,2), 2*bboxDiff(2));
zVec = linspace(bboxPts(1,3),bboxPts(2,3), 2*bboxDiff(3));
[xMat,yMat,zMat] = meshgrid(xVec,yVec,zVec);
normals = faceNormal(myFV.vertices, myFV.faces); % FEX function
in = InPolyedron(myFV.vertices,myFV.faces,normals, [xMat(:),yMat(:),zMat(:)]); % FEX
image3D = reshape(in,size(xMat));
% ALTERNATIVELY, LET'S JUST CREATE SAMPLE INPUT
image3D = false(100,100,100);
image3D(45:55,45:55,:) = true;
image3D(:,45:55,45:55) = true;
image3D(45:55,:,45:55) = true;

% CALCULATE SKELETON
edtImage3D = bwdist(~image3D);
labeledImage = watershed(edtImage3D);
unique(labeledImage(:)) % ans = 1

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <2d36200d-8852-4286-8cc2-4ed4413aec9c@a36g2000yqc.googlegroups.com>...
> Sven:
> Thanks for the test data. But it bombed:
> ??? Undefined function or variable 'xf'.
> Error in ==> test at 16
> 
> Anyway, here is some code to get the medial surface of a 3D volume.
> The algorithm is due to Steve Eddins of the imaging group at the
> Mathworks:
> I start with a 3D image called image3D:
> 
> 		% Start the timer.
> 		[rows cols slices] = size(image3D);
> 		tStart = tic;
> 		% Compute the eroded 3D image because sometimes the particles touch.
> 		SE = strel('ball', 1, 1, 0);
> 		eroded3DImage = imerode(image3D, SE);
> 
> 		% Convert the volume so that the air=1 and the particles=0.
> % 		image3D = 1 - image3D;
> 
> 		% Compute the 3D skeleton using ITK's routine.
> 		% Commented out - didn't seem to work as expected.
> % 		disp('Starting generation of 3D skeleton...');
> % 		image3DSkeleton = ITKSkeleton(image3D);
> 
> 		%------------------------
> 		% Compute the 3D skiz using Steve Eddins proposed method.
> 		% There is a way of using a 3-D watershed transform that I think
> matches closely your inflating balls metaphor.
> 		% Specifically:
> 		% 1.       Compute the distance transform of your 3-D blobs image
> using bwdist.
> 		%          (Fortunately, bwdist runs much faster and takes up less
> memory for the multidimensional case in R2009b than in earlier
> releases.)
> 		% 2.       Compute the watershed transform of the distance transform
> using watershed.
> 		%          The zero-valued pixels in the output of watershed will be
> at the interfaces between the inflated balls.
> 		% I think this is what morphology people call the &#8220;SKIZ,&#8221; or
> &#8220;skeleton by influence zones.&#8221;  The interior of each inflated ball is
> an &#8220;influence zone.&#8221;
> 		% Okay, so here we go......
> 		message = sprintf('Generating 3D distance transform on %d rows by %d
> columns by %d slices 3D image...', rows, cols, slices);
> 		disp(message);
> 		set(handles.txtInfo, 'string', message);
> 		try
> 			%---- Compute the 3D Euclidean Distance Transform.
> 			try
> 				memory
> 				edtImage3D = bwdist(eroded3DImage); % Calculate EDT
> 				clear('eroded3DImage'); % Clear eroded3DImage to free up memory.
> 			catch ME
> 				errorMessage = sprintf('Error running bwdist()\nError Message:\n
> %s', ME.message);
> 				disp(errorMessage);
> 				uiwait(warndlg(errorMessage, 'bwdist error', 'error'));
> 			end
> 			elapsedSeconds = toc(tStart);
> 
> 			%---- Compute the 3D watershed.
> 			message = sprintf('\n%s took %.1f seconds\nGenerating 3D watershed
> transform on %d rows by %d columns by %d slices 3D image...', message,
> elapsedSeconds, rows, cols, slices);
> 			disp(message);
> 			set(handles.txtInfo, 'string', message);
> 			try
> 				memory
> 				labeledImage = watershed_exp(edtImage3D, 6); % Calculate the
> watershed.
> % 				clear('edtImage3D'); % Clear edtImage3D to free up memory.
> 			catch ME
> 				errorMessage = sprintf('Error running watershed()\nError Message:\n
> %s', ME.message);
> 				disp(errorMessage);
> 				uiwait(warndlg(errorMessage, 'watershed error', 'error'));
> 			end
> 			elapsedSeconds = toc(tStart);
> 
> 			%---- Get the pixels with a value of zero
> 			message = sprintf('\n%s %.1f total seconds elapsed so far\nFinding
> zero-value voxels on %d rows by %d columns by %d slices 3D image...',
> message, elapsedSeconds, rows, cols, slices);
> 			disp(message);
> 			image3DSkeleton = (labeledImage == 0);
> 			clear('labeledImage'); % Clear labeledImage to free up memory.
> 
> 
> As you can see, I originally tried a version from ITK but that didn't
> work.  And take note that I called it the "medial surface" and not the
> "skeleton."  Because in 3D the "skeleton" is not necessarily a thin
> single pixel wide line - it can be a surface.  Imagine the space
> between two flat walls.  If that were thinned down, you'd have a sheet
> (surface) lying between the two walls, not a single line.
> 
> Steve's algorithm works well in my application where I wanted the
> medial surface of the void space going between particles that are
> suspended in a 3D volume.
> Good luck,
> -ImageAnalyst
> 

Bahaha... *still* left that typo:

myFV = struct('vertices',xf,'faces',tri); % This is my STL input
should be
myFV = struct('vertices',Xb,'faces',tri); % This is my STL input


"Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <i903o7$8hu$1@fred.mathworks.com>...
> Thanks ImageAnalyst.
> 
> Sorry about the bug. It was a typo that somehow always comes around when you paste in MATLAB code to a browser and then make the tiniest of changes :)  See below which is fixed.
> 
> My question was actually about a triangulated surface, but I am interested in your implementation for a 3D image input. As far as I can tell, the program you gave can be summarised as follows:
> 1. Calculate the bwdist() of your binary 3D image.
> 2. Create a labelled 3D watershed image
> 3. The ridges of the watershed (where label==0) is your skeleton.
> 
> I understand these steps in general, but I found that for the life of me, calling watershed on my distance image simply returns a label matrix full of ones. In your code, you're calling "watershed_exp". Is this a file exchange file, or should the standard IPT "watershed" function work just fine?
> 
> Thanks,
> Sven.
> 
> 
> % CREATE TRIANGULATED INPUT
> load tetmesh
> trep = TriRep(tet, X);
> [tri Xb] = freeBoundary(trep);
> myFV = struct('vertices',xf,'faces',tri); % This is my STL input
> figure, patch(myFV,'EdgeColor','k','FaceColor','r','FaceAlpha',.5), axis image, view(3), camlight 
> 
> % CONVERT TO 3D BINARY IMAGE
> bboxPts = [min(myFV.vertices,[],1);max(myFV.vertices,[],1)];
> bboxDiff = ceil(diff(bboxPts));
> xVec = linspace(bboxPts(1,1),bboxPts(2,1), 2*bboxDiff(1));
> yVec = linspace(bboxPts(1,2),bboxPts(2,2), 2*bboxDiff(2));
> zVec = linspace(bboxPts(1,3),bboxPts(2,3), 2*bboxDiff(3));
> [xMat,yMat,zMat] = meshgrid(xVec,yVec,zVec);
> normals = faceNormal(myFV.vertices, myFV.faces); % FEX function
> in = InPolyedron(myFV.vertices,myFV.faces,normals, [xMat(:),yMat(:),zMat(:)]); % FEX
> image3D = reshape(in,size(xMat));
> % ALTERNATIVELY, LET'S JUST CREATE SAMPLE INPUT
> image3D = false(100,100,100);
> image3D(45:55,45:55,:) = true;
> image3D(:,45:55,45:55) = true;
> image3D(45:55,:,45:55) = true;
> 
> % CALCULATE SKELETON
> edtImage3D = bwdist(~image3D);
> labeledImage = watershed(edtImage3D);
> unique(labeledImage(:)) % ans = 1
> 
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <2d36200d-8852-4286-8cc2-4ed4413aec9c@a36g2000yqc.googlegroups.com>...
> > Sven:
> > Thanks for the test data. But it bombed:
> > ??? Undefined function or variable 'xf'.
> > Error in ==> test at 16
> > 
> > Anyway, here is some code to get the medial surface of a 3D volume.
> > The algorithm is due to Steve Eddins of the imaging group at the
> > Mathworks:
> > I start with a 3D image called image3D:
> > 
> > 		% Start the timer.
> > 		[rows cols slices] = size(image3D);
> > 		tStart = tic;
> > 		% Compute the eroded 3D image because sometimes the particles touch.
> > 		SE = strel('ball', 1, 1, 0);
> > 		eroded3DImage = imerode(image3D, SE);
> > 
> > 		% Convert the volume so that the air=1 and the particles=0.
> > % 		image3D = 1 - image3D;
> > 
> > 		% Compute the 3D skeleton using ITK's routine.
> > 		% Commented out - didn't seem to work as expected.
> > % 		disp('Starting generation of 3D skeleton...');
> > % 		image3DSkeleton = ITKSkeleton(image3D);
> > 
> > 		%------------------------
> > 		% Compute the 3D skiz using Steve Eddins proposed method.
> > 		% There is a way of using a 3-D watershed transform that I think
> > matches closely your inflating balls metaphor.
> > 		% Specifically:
> > 		% 1.       Compute the distance transform of your 3-D blobs image
> > using bwdist.
> > 		%          (Fortunately, bwdist runs much faster and takes up less
> > memory for the multidimensional case in R2009b than in earlier
> > releases.)
> > 		% 2.       Compute the watershed transform of the distance transform
> > using watershed.
> > 		%          The zero-valued pixels in the output of watershed will be
> > at the interfaces between the inflated balls.
> > 		% I think this is what morphology people call the &#8220;SKIZ,&#8221; or
> > &#8220;skeleton by influence zones.&#8221;  The interior of each inflated ball is
> > an &#8220;influence zone.&#8221;
> > 		% Okay, so here we go......
> > 		message = sprintf('Generating 3D distance transform on %d rows by %d
> > columns by %d slices 3D image...', rows, cols, slices);
> > 		disp(message);
> > 		set(handles.txtInfo, 'string', message);
> > 		try
> > 			%---- Compute the 3D Euclidean Distance Transform.
> > 			try
> > 				memory
> > 				edtImage3D = bwdist(eroded3DImage); % Calculate EDT
> > 				clear('eroded3DImage'); % Clear eroded3DImage to free up memory.
> > 			catch ME
> > 				errorMessage = sprintf('Error running bwdist()\nError Message:\n
> > %s', ME.message);
> > 				disp(errorMessage);
> > 				uiwait(warndlg(errorMessage, 'bwdist error', 'error'));
> > 			end
> > 			elapsedSeconds = toc(tStart);
> > 
> > 			%---- Compute the 3D watershed.
> > 			message = sprintf('\n%s took %.1f seconds\nGenerating 3D watershed
> > transform on %d rows by %d columns by %d slices 3D image...', message,
> > elapsedSeconds, rows, cols, slices);
> > 			disp(message);
> > 			set(handles.txtInfo, 'string', message);
> > 			try
> > 				memory
> > 				labeledImage = watershed_exp(edtImage3D, 6); % Calculate the
> > watershed.
> > % 				clear('edtImage3D'); % Clear edtImage3D to free up memory.
> > 			catch ME
> > 				errorMessage = sprintf('Error running watershed()\nError Message:\n
> > %s', ME.message);
> > 				disp(errorMessage);
> > 				uiwait(warndlg(errorMessage, 'watershed error', 'error'));
> > 			end
> > 			elapsedSeconds = toc(tStart);
> > 
> > 			%---- Get the pixels with a value of zero
> > 			message = sprintf('\n%s %.1f total seconds elapsed so far\nFinding
> > zero-value voxels on %d rows by %d columns by %d slices 3D image...',
> > message, elapsedSeconds, rows, cols, slices);
> > 			disp(message);
> > 			image3DSkeleton = (labeledImage == 0);
> > 			clear('labeledImage'); % Clear labeledImage to free up memory.
> > 
> > 
> > As you can see, I originally tried a version from ITK but that didn't
> > work.  And take note that I called it the "medial surface" and not the
> > "skeleton."  Because in 3D the "skeleton" is not necessarily a thin
> > single pixel wide line - it can be a surface.  Imagine the space
> > between two flat walls.  If that were thinned down, you'd have a sheet
> > (surface) lying between the two walls, not a single line.
> > 
> > Steve's algorithm works well in my application where I wanted the
> > medial surface of the void space going between particles that are
> > suspended in a 3D volume.
> > Good luck,
> > -ImageAnalyst
> > 

On Oct 11, 6:40 pm, "Sven" <sven.holco...@gmail.deleteme.com> wrote:
> Thanks ImageAnalyst.
> I understand these steps in general, but I found that for the life of me, calling watershed on my distance image simply returns a label matrix full of ones. In your code, you're calling "watershed_exp". Is this a file exchange file, or should the standard IPT "watershed" function work just fine?
---------------------------------------------------------------------------
Sorry about that.  watershed_exp() is a new experimental version of watershed that Steve was working on.  I'm not sure of the status of it, like if it made it into R2010b or not.  You'd have to contact him if you want it, as I wouldn't feel authorized to give it to you.