Skip to content

A Pictorial Guide to TEASAR Skeletonization

Jump to bottom
William Silversmith edited this page Aug 8, 2019 · 18 revisions

Introduction

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.

Short Narrative

This algorithm works by finding a root point on a 3D object and then serially tracing paths via dijksta's shortest path algorithm through a penalty field to the most distant unvisited point. After each pass, there is a sphere (really implemented as a circumscribing cube) that expands around each vertex in the current path that marks part of the object as visited.

Guide

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:

Euclidean distance transform of an isotropic dimension square. The square is colored in a gradient that is dark at the edges and bright in the center. There is a fringe of brightness along the diagonals.
A dumbbell shape colored with stippling to show a darker rim and a lighter interior.

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

Showing PDRF equation.

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 (and a deviation from TEASAR). The typical value for k is about 100,000. The typical range for r is between 1 and 16 and is defaulted to 4.

Showing augmented PDRF equation.

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). This is a deviation from TEASAR.

Dr is a field of geodesic distances from the root.

Drawing a path from point s (the root) to point t through the original Y-fork shape.

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. In original TEASAR, the selection of the target involved the PDRF, but we found the results were more sensible by using just the distance from source field.

The region of the shape (and some outside of the shape) surrounding the previously drawn path is invalidated.

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.

A path is drawn from root s to a new point t on the other branch of the Y-fork that was not invalidated.

7. Draw a path to the next surviving max distance point and invalidate. Repeat until no foreground voxels survive.

Clone this wiki locally