🚀 go-pugleaf

Well, I've only gotten a couple of nibbles for 3D interrogation
software, but I have had a bunch of folk ask me about my 3D medial
axis algorithm.  So, here it is.  I can't believe this isn't out
there somewhere in the literature -- or else it doesn't work
nearly as well as I think it does.  Now I'm paranoid.
Anyway, take a look. No guarantees.  

My method is pretty brute force, and works only for binary images,
is slow, but seems to work OK, given lots of time and effort.  It is an
extension of something John Gauch did for 2D images when he was
a grad student at UNC.
  
It consists of:

1) Find the 3D radius functions for all voxels (this is done simply
by an erosive method).

2) Sort the non-zero voxels by radius.

3) For each non-zero voxel, starting with those with the larger radii{
	For each voxel with a smaller radius (e.g. the rest of the list){
		if the voxel with the larger radius "covers" the voxel
			with the smaller radius, discard the smaller one.
		}
	}

4) What's left is on the axis.


Covering is simple to calculate, e.g. let X be the center of the
"big" voxel with radius R and x be the center of the "little" voxel with
radius r, then the big one covers the little one if the distance
|X-x|+r is less than or equal to R.

This leaves some bad spots due to discretization, where there is a loss 
of connectivity at some places where the radius decreases.  One example is
a configuration like :

*****
**********
*****

which has a radius function of:

11111
2222111111
11111

(the rest of the object is off to the left) and which gives a skeleton of:

**** *****

instead of:  

**********

.

Thus, I have a few heuristics to go back and look to fill in missing
spots.  One of the reasons that I am looking for a 3D interrogation 
tool is so I can see exactly how well this is working on complicated 
data sets.

billo