Kimimaro
Skeletonize densely labeled 3D image segmentations with TEASAR. (Medial Axis Transform)
Install / Use
/learn @seung-lab/KimimaroREADME
Kimimaro: Skeletonize Densely Labeled Images
# Produce SWC files from volumetric images.
kimimaro forge labels.npy --progress # writes to ./kimimaro_out/
kimimaro view kimimaro_out/10.swc
Rapidly skeletonize all non-zero labels in 2D and 3D numpy arrays using a TEASAR derived method. The returned list of skeletons is in the format used by cloud-volume. A skeleton is a stick figure 1D representation of a 2D or 3D object that consists of a graph of verticies linked by edges. A skeleton where the verticies also carry a distance to the nearest boundary they were extracted from is called a "Medial Axis Transform", which Kimimaro provides.
Skeletons are a compact representation that can be used to visualize objects, trace the connectivity of an object, or otherwise analyze the object's geometry. Kimimaro was designed for use with high resolution neurons extracted from electron microscopy data via AI segmentation, but it can be applied to many different fields.
On an Apple Silicon M1 arm64 chip (Firestorm cores 3.2 GHz max frequency), this package processed a 512x512x100 volume with 333 labels in 20 seconds. It processed a 512x512x512 volume (connectomics.npy) with 2124 labels in 187 seconds.
pip Installation
If a binary is available for your platform:
pip install kimimaro
# installs additional libraries to accelerate some
# operations like join_close_components
pip install "kimimaro[accel]"
# Makes the kimimaro view command work
pip install "kimimaro[view]"
# Enables TIFF generation on the CLI
pip install "kimimaro[tif]"
# Enables reading NIBABEL, NRRD, TIFF, CRACKLE on the CLI
pip install "kimimaro[all_formats]"
# Install all optional dependencies
pip install "kimimaro[all]"
Otherwise, you'll also need a C++ compiler:
sudo apt-get install python3-dev g++ # ubuntu linux
Example
<p style="font-style: italics;" align="center"> <img height=512 src="https://raw.githubusercontent.com/seung-lab/kimimaro/master/kimimaro_512x512x512_benchmark.png" alt="A Densely Labeled Volume Skeletonized with Kimimaro" /><br> Fig. 2: Memory Usage on a 512x512x512 Densely Labeled Volume (`connectomics.npy`) </p>Figure 2 shows the memory usage and processessing time (~390 seconds, about 6.5 minutes) required when Kimimaro 1.4.0 was applied to a 512x512x512 cutout, labels, from a connectomics dataset, connectomics.npy containing 2124 connected components. The different sections of the algorithm are depicted. Grossly, the preamble runs for about half a minute, skeletonization for about six minutes, and finalization within seconds. The peak memory usage was about 4.5 GB. The code below was used to process labels. The processing of the glia was truncated in due to a combination of fix_borders and max_paths.
Kimimaro has come a long way. Version 0.2.1 took over 15 minutes and had a Preamble run time twice as long on the same dataset.
On a Macbook Pro M3, the same settings now complete in 94 seconds (1.6 minutes) on version 5.4.0. With xs3d 1.11.0, cross section analysis takes 215 seconds (3.6 minutes).
Python Interface
# LISTING 1: Producing Skeletons from a labeled image.
import kimimaro
# To obtain this 512 MB segmentation sample volume:
# pip install crackle-codec
import crackle
labels = crackle.load("benchmarks/connectomics.npy.ckl.gz")
skels = kimimaro.skeletonize(
labels,
teasar_params={
"scale": 1.5,
"const": 300, # physical units
"pdrf_scale": 100000,
"pdrf_exponent": 4,
"soma_acceptance_threshold": 3500, # physical units
"soma_detection_threshold": 750, # physical units
"soma_invalidation_const": 300, # physical units
"soma_invalidation_scale": 2,
"max_paths": 300, # default None
},
# object_ids=[ ... ], # process only the specified labels
# extra_targets_before=[ (27,33,100), (44,45,46) ], # target points in voxels
# extra_targets_after=[ (27,33,100), (44,45,46) ], # target points in voxels
dust_threshold=1000, # skip connected components with fewer than this many voxels
anisotropy=(16,16,40), # default True
fix_branching=True, # default True
fix_borders=True, # default True
fill_holes=False, # default False
fix_avocados=False, # default False
progress=True, # default False, show progress bar
parallel=1, # <= 0 all cpu, 1 single process, 2+ multiprocess
parallel_chunk_size=100, # how many skeletons to process before updating progress bar
)
# LISTING 2: Combining skeletons produced from
# adjacent or overlapping images.
import kimimaro
from osteoid import Skeleton
skels = ... # a set of skeletons produced from the same label id
skel = Skeleton.simple_merge(skels).consolidate()
skel = kimimaro.postprocess(
skel,
dust_threshold=1000, # physical units
tick_threshold=3500 # physical units
)
# LISTING 3: Adding cross sectional area to skeletons
# Cross section planes are defined by normal vectors. Those
# vectors come from the difference between adjacent vertices.
skels = ... # one or more skeletons produced from a single image
skels = kimimaro.cross_sectional_area(
labels, skels,
anisotropy=(16,16,40),
smoothing_window=5, # rolling average window of plane normals
progress=True,
)
skel = skels[0]
skel.cross_sectional_area # array of cross sectional areas
skel.cross_sectional_area_contacts # non-zero contacted the image border
# Split input skeletons into connected components and
# then join the two nearest vertices within `radius` distance
# of each other until there is only a single connected component
# or no pairs of points nearer than `radius` exist.
# Fuse all remaining components into a single skeleton.
skel = kimimaro.join_close_components([skel1, skel2], radius=1500) # 1500 units threshold
skel = kimimaro.join_close_components([skel1, skel2], radius=None) # no threshold
# Given synapse centroids (in voxels) and the SWC integer label you'd
# like to assign (e.g. for pre-synaptic and post-synaptic) this finds the
# nearest voxel to the centroid for that label.
# Input: { label: [ ((x,y,z), swc_label), ... ] }
# Returns: { (x,y,z): swc_label, ... }
extra_targets = kimimaro.synapses_to_targets(labels, synapses)
# LISTING 4: Drawing a centerline between
# preselected points on a binary image.
# This is a much simpler option for when
# you know exactly what you want, but may
# be less efficient for large scale procesing.
skel = kimimaro.connect_points(
labels == 67301298,
start=(3, 215, 202),
end=(121, 426, 227),
anisotropy=(32,32,40),
)
# LISTING 5: Using skeletons to oversegment existing
# segmentations for integration into proofreading systems
# that on merging atomic labels. oversegmented_labels
# is returned numbered from 1. skels is a copy returned
# with the property skel.segments that associates a label
# to each vertex (labels will not be unique if downsampling
# is used)
oversegmented_labels, skels = kimimaro.oversegment(
labels, skels,
anisotropy=(32,32,40),
downsample=10,
)
connectomics.npy is multilabel connectomics data derived from pinky40, a 2018 experimental automated segmentation of ~1.5 million cubic micrometers of mouse visual cortex. It is an early predecessor to the now public pinky100_v185 segmentation that can be found at https://microns-explorer.org/phase1 You will need to run lzma -d connectomics.npy.lzma to obtain the 512x512x512 uint32 volume at 32x32x40 nm<sup>3</sup> resolution.
CLI Interface
The CLI supports producing skeletons from a single image as SWCs and viewing the resulting SWC files one at a time. By default, the SWC files are written to ./kimimaro_out/$LABEL.swc.
Here's an equivalent example to the code above.
kimimaro forge labels.npy --scale 4 --const 10 --soma-detect 1100 --soma-accept 3500 --soma-scale 1 --soma-const 300 --anisotropy 16,16,40 --fix-borders --progress
Visualize the your data:
kimimaro view 1241241.swc # visualize skeleton
kimimaro view labels.npy # visualize segmentation
It can also convert binary image skeletons produced by thinning algorithms into SWC files and back. This can be helpful for comparing different skeletonization algorithms or even just using their results.
kimimaro swc from binary_image.tiff # -> binary_image.swc
kimimaro swc to --format tiff binary_image.swc # -> binary_image.tiff or npy
Tweaking kimimaro.skeletonize Parameters
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 a circumscribing cube) that expands around each vertex in the current path that marks part of the object as visited.
For a visual tutorial on the basics of the skeletonization procedure, check out this wiki article: A Pictorial Guide to TEASAR Skeletonization
For more detailed information, read below or the TEASAR paper (though we deviate from TEASAR in a few places). [1]
scale and const
Usually, the most important parameters to tweak are scale and const which control the radius of this invalidation sphere according to the equation `r(x,y,z) = scale * DBF(x,y,z
