-
Notifications
You must be signed in to change notification settings - Fork 25
A Pictorial Guide to TEASAR Skeletonization
Below is an illustrated guide for the basic outline of how Kimimaro skeletonizes a label in a 3D image. For clarity, this guide does not necessarily present the steps in the same order as the code does. There are some deviations from base TEASAR which will be noted. We also do not cover the important effects of e.g. fix_branches
or fix_borders
, nor do we cover the Preamble or handling of multiple labels. The idea is to give you a basic intuition for the algorithm so that the more detailed descriptions become intelligible.
If this guide is helpful for people, I might expand it.
1. TEASAR Skeletonization 2D Examples. You have a shape. A shape with a Y-fork in it is depicted.
2. Raster scan to find a foreground voxel. This finds an arbitrary foreground voxel.
3. Find the furthest point from this starting voxel. This is the root.
4a. Create a penalty field to guide skeletons through the center of the shape. D(x,y,z) is the euclidean distance of each pixel from the boundary. (L2 norm).
i.e. we are computing the Euclidean Distance Transform. For example, here is the distance transform of a square:
4b. Create a penalty field to guide skeletons through the center of the shape. P(x,y,z) is a field that is valued low in the center of an object and high at the edges. (Darker is higher value).
4c. Showing the equation of the Penalized Distance from Root Field (PDRF). `P(x,y,z) = k * (1 - D / max(D)^1.01) ^ r` where k and r are the user specified constants `pdrf_scale` and `pdrf_exp` respectively. 1.01 is arbitrary. The typical value for k is about 100,000. The typical range for r is between 1 and 16 and is defaulted to 4.
Sometimes the value of P(x,y,z) can collapse lower than float32 precision which would cause the skeleton path to wander randomly. We add a small gradient to provide direction in those cases. P(x,y,z) = k * (1 - D / max(D)^1.01) ^ r + Dr / max(Dr)
Dr is a field of geodesic distances from the root.
5. Use Dijkstra's shortest path algorithm to draw a path from the root (s) to the target (t), which is the most geodesicaly distant point from s, through the penalty field P.
6. Invalidate foregound voxels by rolling a ball of dynamic radius r(x,y,z)
down the path.
r(x,y,z) = scale * D(x,y,z) + const
where D(x,y,z)
is the value of the distance transform at that voxel. scale
and const
are user provided.
7. Draw a path to the next surviving max distance point and invalidate. Repeat until no foreground voxels survive.