{
"info": {
"author": "William Silversmith, Alex Bae, Forrest Collman",
"author_email": "ws9@princeton.edu",
"bugtrack_url": null,
"classifiers": [
"Development Status :: 4 - Beta",
"Intended Audience :: Developers",
"License :: OSI Approved :: GNU General Public License v3 (GPLv3)",
"Programming Language :: Python :: 2",
"Programming Language :: Python :: 2.7",
"Programming Language :: Python :: 3",
"Programming Language :: Python :: 3.5",
"Programming Language :: Python :: 3.6",
"Programming Language :: Python :: 3.7",
"Topic :: Scientific/Engineering"
],
"description": "[](https://travis-ci.org/seung-lab/kimimaro) [](https://badge.fury.io/py/kimimaro) \n\n# Kimimaro: Skeletonize Densely Labeled Images\n\nRapidly 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](https://github.com/seung-lab/cloud-volume/wiki/Advanced-Topic:-Skeletons). \n\nOn a 3.7 GHz Intel i7 processor, this package processed a 512x512x100 volume with 333 labels in under a minute. It processed a 512x512x512 volume with 2124 labels in eight to thirteen minutes (depending on whether `fix_branching` is set).\n\n
\n \nFig. 1: A Densely Labeled Volume Skeletonized with Kimimaro\n
\n\n## `pip` Installation \n\n*Requires C++ compiler.*\n\n```bash\nsudo apt-get install python3-dev g++\npip3 install numpy\npip3 install kimimaro \n```\n\nIn the future, we may create a fully binary distribution. \n\n## Example\n\n
\n \nFig. 2: Memory Usage on a 512x512x512 Densely Labeled Volume\n
\n\nFigure 2 shows the memory usage and processessing time (~390 seconds, about 6.5 minutes) required when Kimimaro 0.5.2 was applied to a 512x512x512 cutout, *labels*, from a connectomics dataset 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.1 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. \n\nKimimaro 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. \n\n```python\n# LISTING 1: Producing Skeletons from a labeled image.\n\nimport kimimaro\n\nlabels = np.load(...)\n\nskels = kimimaro.skeletonize(\n labels, \n teasar_params={\n 'scale': 4,\n 'const': 500, # physical units\n 'pdrf_exponent': 4,\n 'pdrf_scale': 100000,\n 'soma_detection_threshold': 1100, # physical units\n 'soma_acceptance_threshold': 3500, # physical units\n 'soma_invalidation_scale': 1.0,\n 'soma_invalidation_const': 300, # physical units\n 'max_paths': 50, # default None\n },\n # object_ids=[ ... ], # process only the specified labels\n # extra_targets_before=[ (27,33,100), (44,45,46) ], # target points in voxels\n # extra_targets_after=[ (27,33,100), (44,45,46) ], # target points in voxels\n dust_threshold=1000, # skip connected components with fewer than this many voxels\n anisotropy=(16,16,40), # default True\n fix_branching=True, # default True\n fix_borders=True, # default True\n progress=True, # default False, show progress bar\n parallel=1, # <= 0 all cpu, 1 single process, 2+ multiprocess\n parallel_chunk_size=100, # how many skeletons to process before updating progress bar\n)\n\n# LISTING 2: Combining skeletons produced from \n# adjacent or overlapping images.\n\nimport kimimaro\nfrom cloudvolume import PrecomputedSkeleton\n\nskels = ... # a set of skeletons produced from the same label id\nskel = PrecomputedSkeleton.simple_merge(skels).consolidate()\nskel = kimimaro.postprocess(\n skel, \n dust_threshold=1000, # physical units\n tick_threshold=3500 # physical units\n)\n\n# Split input skeletons into connected components and\n# then join the two nearest vertices within `radius` distance\n# of each other until there is only a single connected component\n# or no pairs of points nearer than `radius` exist. \n# Fuse all remaining components into a single skeleton.\nskel = kimimaro.join_close_components([skel1, skel2], radius=1500) # 1500 units threshold\nskel = kimimaro.join_close_components([skel1, skel2], radius=None) # no threshold\n```\n\n### Tweaking `kimimaro.skeletonize` Parameters\n\nThis 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. \n\nFor a visual tutorial on the basics of the skeletonization procedure, check out this wiki article: [A Pictorial Guide to TEASAR Skeletonization](https://github.com/seung-lab/kimimaro/wiki/A-Pictorial-Guide-to-TEASAR-Skeletonization)\n\nFor more detailed information, [read below](https://github.com/seung-lab/kimimaro#ii-skeletonization) or the [TEASAR paper](https://ieeexplore.ieee.org/abstract/document/883951/) (though we [deviate from TEASAR](https://github.com/seung-lab/kimimaro#teasar-derived-algorthm) in a few places). [1]\n\n#### `scale` and `const`\n\nUsually, 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) + const` where the dimensions are physical (e.g. nanometers, i.e. corrected for anisotropy). `DBF(x,y,z)` is the physical distance from the shape boundary at that point. \n\nCheck out this [wiki article](https://github.com/seung-lab/kimimaro/wiki/Intuition-for-Setting-Parameters-const-and-scale) to help refine your intuition.\n\n#### `anisotropy`\n\nRepresents the physical dimension of each voxel. For example, a connectomics dataset might be scanned with an electron microscope at 4nm x 4nm per pixel and stacked in slices 40nm thick. i.e. `anisotropy=(4,4,40)`. You can use any units so long as you are consistent.\n\n#### `dust_threshold`\n\nThis threshold culls connected components that are smaller than this many voxels. \n\n#### `extra_targets_after` and `extra_targets_before` \n\n`extra_targets_after` provides additional voxel targets to trace to after the morphological tracing algorithm completes. For example, you might add known synapse locations to the skeleton. \n\n`extra_targets_before` is the same as `extra_targets_after` except that the additional targets are front-loaded and the paths that they cover are invalidated. This may affect the results of subsequent morphological tracing.\n\n#### `max_paths` \n\nLimits the number of paths that can be drawn for the given label. Certain cells, such as glia, that may not be important for the current analysis may be expensive to process and can be aborted early. \n\n#### `pdrf_scale` and `pdrf_exponent`\n\nThe `pdrf_scale` and `pdrf_exponent` represent parameters to the penalty equation that takes the euclidean distance field (**D**) and augments it so that cutting closer to the border is very penalized to make dijkstra take paths that are more centered. \n\nPr = `pdrf_scale` * (1 - **D** / max(**D**)) `pdrf_exponent` + (directional gradient < 1.0). \n\nThe default settings should work fairly well, but under large anisotropies or with cavernous morphologies, it's possible that you might need to tweak it. If you see the skeleton go haywire inside a large area, it could be a collapse of floating point precision. \n\n#### `soma_acceptance_threshold` and `soma_detection_threshold`\n\nWe process somas specially because they do not have a tubular geometry and instead should be represented in a hub and spoke manner. `soma_acceptance_threshold` is the physical radius (e.g. in nanometers) beyond which we classify a connected component of the image as containing a soma. The distance transform's output is depressed by holes in the label, which are frequently produced by segmentation algorithms on somata. We can fill them, but the hole filling algorithm we use is slow so we would like to only apply it occasionally. Therefore, we set a lower threshold, the `soma_acceptance_threshold`, beyond which we fill the holes and retest the soma. \n\n#### `soma_invalidation_scale` and `soma_invalidation_const` \n\nOnce we have classified a region as a soma, we fix root of the skeletonization algorithm at one of the points of maximum distance from the boundary (usually there is only one). We then mark as visited all voxels around that point in a spherical radius described by `r(x,y,z) = soma_invalidation_scale * DBF(x,y,z) + soma_invalidation_const` where DBF(x,y,z) is the physical distance from the shape boundary at that point. If done correctly, this can prevent skeletons from being drawn to the boundaries of the soma, and instead pulls the skeletons mainly into the processes extending from the cell body. \n\n#### `fix_borders`\n\nThis feature makes it easier to connect the skeletons of adjacent image volumes that do not fit in RAM. If enabled, skeletons will be deterministically drawn to the approximate center of the 2D contact area of each place where the shape contacts the border. This can affect the performance of the operation positively or negatively depending on the shape and number of contacts. \n\n#### `fix_branching` \n\nYou'll probably never want to disable this, but base TEASAR is infamous for forking the skeleton at branch points way too early. This option makes it preferential to fork at a more reasonable place at a significant performance penalty. \n\n#### `progress`\n\nShow a progress bar once the skeletonization phase begins.\n\n#### `parallel` \n\nUse a pool of processors to skeletonize faster. Each process allocatable task is the skeletonization of one connected component (so it won't help with a single label that takes a long time to skeletonize). This option also affects the speed of the initial euclidean distance transform, which is parallel enabled and is the most expensive part of the Preamble (described below). \n\n#### `parallel_chunk_size` \n\nThis only applies when using parallel. This sets the number of skeletons a subprocess will extract before returning control to the main thread, updating the progress bar, and acquiring a new task. If this value is set too low (e.g. < 10-20) the cost of interprocess communication can become significant and even dominant. If it is set too high, task starvation may occur for the other subprocesses if a subprocess gets a particularly hard skeleton and they complete quickly. Progress bar updates will be infrequent if the value is too high as well. \n\nThe actual chunk size used will be `min(parallel_chunk_size, len(cc_labels) // parallel)`. `cc_labels` represents the number of connected components in the sample. \n\n### Performance Tips\n\n- If you only need a few labels skeletonized, pass in `object_ids` to bypass processing all the others. If `object_ids` contains only a single label, the masking operation will run faster.\n- You may save on peak memory usage by using a `cc_safety_factor` < 1, only if you are sure the connected components algorithm will generate many fewer labels than there are pixels in your image.\n- Larger TEASAR parameters scale and const require processing larger invalidation regions per path.\n- Set `pdrf_exponent` to a small power of two (e.g. 1, 2, 4, 8, 16) for a small speedup.\n- If you are willing to sacrifice the improved branching behavior, you can set `fix_branching=False` for a moderate 1.1x to 1.5x speedup (assuming your TEASAR parameters and data allow branching).\n- If your dataset contains important cells (that may in fact be the seat of consciousness) but they take significant processing power to analyze, you can save them to savor for later by setting `max_paths` to some reasonable level which will abort and proceed to the next label after the algorithm detects that that at least that many paths will be needed.\n- Parallel distributes work across connected components and is generally a good idea if you have the cores and memory. Not only does it make single runs proceed faster, but you can also practically use a much larger context; that improves soma processing as they are less likely to be cut off. The Preamble of the algorithm (detailed below) is still single threaded at the moment, so task latency increases with size. \n- If `parallel_chunk_size` is set very low (e.g. < 10) during parallel operation, interprocess communication can become a significant overhead. Try raising this value. \n\n## Motivation\n\nThe connectomics field commonly generates very large densely labeled volumes of neural tissue. Skeletons are one dimensional representations of two or three dimensional objects. They have many uses, a few of which are visualization of neurons, calculating global topological features, rapidly measuring electrical distances between objects, and imposing tree structures on neurons (useful for computation and user interfaces). There are several ways to compute skeletons and a few ways to define them [4]. After some experimentation, we found that the TEASAR [1] approach gave fairly good results. Other approaches include topological thinning (\"onion peeling\") and finding the centerline described by maximally inscribed spheres. Ignacio Arganda-Carreras, an alumnus of the Seung Lab, wrote a topological thinning plugin for Fiji called [Skeletonize3d](https://imagej.net/Skeletonize3D). \n\nThere are several implementations of TEASAR used in the connectomics field [3][5], however it is commonly understood that implementations of TEASAR are slow and can use tens of gigabytes of memory. Our goal to skeletonize all labels in a petavoxel scale image quickly showed clear that existing sparse implementations are impractical. While adapting a sparse approach to a cloud pipeline, we noticed that there are inefficiencies in repeated evaluation of the Euclidean Distance Transform (EDT), the repeated evaluation of the connected components algorithm, in the construction of the graph used by Dijkstra's algorithm where the edges are implied by the spatial relationships between voxels, in the memory cost, quadratic in the number of voxels, of representing a graph that is implicit in image, in the unnecessarily large data type used to represent relatively small cutouts, and in the repeated downloading of overlapping regions. We also found that the naive implmentation of TEASAR's \"rolling invalidation ball\" unnecessarily reevaluated large numbers of voxels in a way that could be loosely characterized as quadratic in the skeleton path length. \n\nWe further found that commodity implementations of the EDT supported only binary images. We were unable to find any available Python or C++ libraries for performing Dijkstra's shortest path on an image. Commodity implementations of connected components algorithms for images supported only binary images. Therefore, several libraries were devised to remedy these deficits (see Related Projects). \n\n## Why TEASAR?\n\nTEASAR: Tree-structure Extraction Algorithm for Accurate and Robust skeletons, a 2000 paper by M. Sato and others [1], is a member of a family of algorithms that transform two and three dimensional structures into a one dimensional \"skeleton\" embedded in that higher dimension. One might concieve of a skeleton as extracting a stick figure drawing from a binary image. This problem is more difficult than it might seem. There are different situations one must consider when making such a drawing. For example, a stick drawing of a banana might merely be a curved centerline and a drawing of a doughnut might be a closed loop. In our case of analyzing neurons, sometimes we want the skeleton to include spines, short protrusions from dendrites that usually have synapses attached, and sometimes we want only the characterize the run length of the main trunk of a neurite. \n\nAdditionally, data quality issues can be challenging as well. If one is skeletonizing a 2D image of a doughnut, but the angle were sufficiently declinated from the ring's orthogonal axis, would it even be possible to perform this task accurately? In a 3D case, if there are breaks or mergers in the labeling of a neuron, will the algorithm function sensibly? These issues are common in both manual and automatic image sementations.\n\nIn our problem domain of skeletonizing neurons from anisotropic voxel labels, our chosen algorithm should produce tree structures, handle fine or coarse detail extraction depending on the circumstances, handle voxel anisotropy, and be reasonably efficient in CPU and memory usage. TEASAR fufills these criteria. Notably, TEASAR doesn't guarantee the centeredness of the skeleton within the shape, but it makes an effort. The basic TEASAR algorithm is known to cut corners around turns and branch too early. A 2001 paper by members of the original TEASAR team describes a method for reducing the early branching issue on page 204, section 4.2.2. [2]\n\n## TEASAR Derived Algorthm\n\nWe implemented TEASAR but made several deviations from the published algorithm in order to improve path centeredness, increase performance, handle bulging cell somas, and enable efficient chunked evaluation of large images. We opted not to implement the gradient vector field step from [2] as our implementation is already quite fast. The paper claims a reduction of 70-85% in input voxels, so it might be worth investigating. \n\nIn order to work with images that contain many labels, our general strategy is to perform as many actions as possible in such a way that all labels are treated in a single pass. Several of the component algorithms (e.g. connected components, euclidean distance transform) in our implementation can take several seconds per a pass, so it is important that they not be run hundreds or thousands of times. A large part of the engineering contribution of this package lies in the efficiency of these operations which reduce the runtime from the scale of hours to minutes. \n\nGiven a 3D labeled voxel array, *I*, with N >= 0 labels, and ordered triple describing voxel anisotropy *A*, our algorithm can be divided into three phases, the pramble, skeletonization, and finalization in that order.\n\n### I. Preamble\n\nThe Preamble takes a 3D image containing *N* labels and efficiently generates the connected components, distance transform, and bounding boxes needed by the skeletonization phase.\n\n1. To enhance performance, if *N* is 0 return an empty set of skeletons.\n2. Label the M connected components, *Icc*, of *I*.\n3. To save memory, renumber the connected components in order from 1 to *M*. Adjust the data type of the new image to the smallest uint type that will contain *M* and overwrite *Icc*.\n4. Generate a mapping of the renumbered *Icc* to *I* to assign meaningful labels to skeletons later on and delete *I* to save memory.\n5. Compute *E*, the multi-label anisotropic Euclidean Distance Transform of *Icc* given *A*. *E* treats all interlabel edges as transform edges, but not the boundaries of the image. Black pixels are considered background.\n6. Gather a list, *Lcc* of unique labels from *Icc* and threshold which ones to process based on the number of voxels they represent to remove \"dust\".\n7. In one pass, compute the list of bounding boxes, *B*, corresponding to each label in *Lcc*.\n\n### II. Skeletonization \n\nIn this phase, we extract the tree structured skeleton from each connected component label. Below, we reference variables defined in the Preamble. For clarity, we omit the soma specific processing and hold `fix_branching=True`. \n\nFor each label *l* in *Lcc* and *B*...\n\n1. Extract *Il*, the cropped binary image tightly enclosing *l* from *Icc* using *Bl*\n2. Using *Il* and *Bl*, extract *El* from *E*. *El* is the cropped tightly enclosed EDT of *l*. This is much faster than recomputing the EDT for each binary image.\n3. Find an arbitrary foreground voxel and using that point as a source, compute the anisotropic euclidean distance field for *Il*. The coordinate of the maximum value is now \"the root\" *r*.\n4. From *r*, compute the euclidean distance field and save it as the distance from root field *Dr*.\n5. Compute the penalized distance from root field *Pr* = `pdrf_scale` * ((1 - *El* / max(*El*)) ^ `pdrf_exponent`) + *Dr* / max(*Dr*). \n6. While *Il* contains foreground voxels:\n 1. Identify a target coordinate, *t*, as the foreground voxel with maximum distance in *Dr* from *r*.\n 2. Draw the shortest path *p* from *r* to *t* considering the voxel values in *Pr* as edge weights.\n 3. For each vertex *v* in *p*, extend an invalidation cube of physical side length computed as `scale` * *El*(*v*) + `const` and convert any foreground pixels in *Il* that overlap with these cubes to background pixels.\n 4. (Only if `fix_branching=True`) For each vertex coordinate *v* in *p*, set *Pr*(*v*) = 0.\n 5. Append *p* to a list of paths for this label.\n7. Using *El*, extract the distance to the nearest boundary each vertex in the skeleton represents.\n8. For each raw skeleton extracted from *Il*, translate the vertices by *Bl* to correct for the translation the cropping operation induced.\n9. Multiply the vertices by the anisotropy *A* to place them in physical space.\n\nIf soma processing is considered, we modify the root (*r*) search process as follows: \n\n1. If max(*El*) > `soma_detection_threshold`...\n 1. Fill toplogical holes in *Il*. Soma are large regions that often have dust from imperfect automatic labeling methods.\n 2. Recompute *El* from this cleaned up image.\n 3. If max(*El*) > `soma_acceptance_threshold`, divert to soma processing mode.\n2. If in soma processing mode, continue, else go to step 3 in the algorithm above.\n3. Set *r* to the coordinate corresponding to max(*El*)\n4. Create an invalidation sphere of physical radius `soma_invalidation_scale` * max(*El*) + `soma_invalidation_const` and erase foreground voxels from *Il* contained within it. This helps prevent errant paths from being drawn all over the soma.\n5. Continue from step 4 in the above algorithm.\n\n### III. Finalization\n\nIn the final phase, we agglomerate the disparate connected component skeletons into single skeletons and assign their labels corresponding to the input image. This step is artificially broken out compared to how intermingled its implementation is with skeletonization, but it's conceptually separate.\n\n## Deviations from TEASAR\n\nThere were several places where we took a different approach than called for by the TEASAR authors.\n\n### Using DAF for Targets, PDRF for Pathfinding\n\nThe original TEASAR algorithm defines the Penalized Distance from Root voxel Field (PDRF, *Pr* above) as:\n\n```\nPDRF = 5000 * (1 - DBF / max(DBF))^16 + DAF\n```\n\nDBF is the Distance from Boundary Field (*El* above) and DAF is the Distance from Any voxel Field (*Dr* above). \n\nWe found the addition of the DAF tended to perturb the skeleton path from the centerline better described by the inverted DBF alone. We also found it helpful to modify the constant and exponent to tune cornering behavior. Initially, we completely stripped out the addition of the DAF from the PDRF, but this introduced a different kind of problem. The exponentiation of the PDRF caused floating point values to collapse in wide open spaces. This made the skeletons go crazy as they traced out a path described by floating point errors. \n\nThe DAF provides a very helpful gradient to follow between the root and the target voxel, we just don't want that gradient to knock the path off the centerline. Therefore, in light of the fact that the PDRF base field is very large, we add the normalized DAF which is just enough to overwhelm floating point errors and provide direction in wide tubes and bulges. \n\nThe original paper also called for selecting targets using the max(PDRF) foreground values. However, this is a bit strange since the PDRF values are dominated by boundary effects rather than a pure distance metric. Therefore, we select targets from the max(DAF) forground value.\n\n### Zero Weighting Previous Paths (`fix_branching=True`)\n\nThe 2001 skeletonization paper [2] called for correcting early forking by computing a DAF using already computed path vertices as field sources. This allows Dijkstra's algorithm to trace the existing path cost free and diverge from it at a closer point to the target. \n\nAs we have strongly deemphasized the role of the DAF in dijkstra path finding, computing this field is unnecessary and we only need to set the PDRF to zero along the path of existing skeletons to achieve this effect. This saves us an expensive repeated DAF calculation per path. \n\nHowever, we still incur a substantial cost for taking this approach because we had been computing a dijkstra \"parental field\" that recorded the shortest path to the root from every foreground voxel. We then used this saved result to rapidly compute all paths. However, as this zero weighting modification makes successive calculations dependent upon previous ones, we need to compute Dijkstra's algorithm anew for each path.\n\n### Non-Overlapped Chunked Processing (`fix_borders=True`)\n\nWhen processing large volumes, a sensible approach for mass producing skeletons is to chunk the volume, process the chunks independently, and merge the resulting skeleton fragments at the end. However, this is complicated by the \"edge effect\" induced by a loss of context which makes it impossible to expect the endpoints of skeleton fragments produced by adjacent chunks to align. In contrast, it is easy to join mesh fragments because the vertices of the edge of mesh fragments lie at predictable identical locations given one pixel of overlap. \n\nPreviously, we had used 50% overlap to join adjacent skeleton fragments which increased the compute cost of skeletonizing a large volume by eight times. However, if we could force skeletons to lie at predictable locations on the border, we could use single pixel overlap and copy the simple mesh joining approach. As an (incorrect but useful) intuition for how one might go about this, consider computing the centroid of each connected component on each border plane and adding that as a required path target. This would guarantee that both sides of the plane connect at the same pixel. However, the centroid may not lie inside of non-convex hulls so we have to be more sophisticated and select some real point inside of the shape.\n\nTo this end, we again repurpose the euclidean distance transform and apply it to each of the six planes of connected components and select the maximum value as a mandatory target. This works well for many types of objects that contact a single plane and have a single maximum. However, we must treat the corners of the box and shapes that have multiple maxima. \n\nTo handle shapes that contact multiple sides of the box, we simply assign targets to all connected components. If this introduces a cycle in post-processing, we already have cycle removing code to handle it in Igneous. If it introduces tiny useless appendages, we also have code to handle this. \n\nIf a shape has multiple distance transform maxima, it is important to choose the same pixel without needing to communicate between spatially adjacent tasks which may run at different times on different machines. Additionally, the same plane on adjacent tasks has the coordinate system flipped. One simple approach might be to pick the coordinate with minimum x and y (or some other coordinate based criterion) in one of the coordinate frames, but this requires tracking the flips on all six planes and is annoying. Instead, we use a series of coordinate-free topology based filters which is both more fun, effort efficient, and picks something reasonable looking. A valid criticism of this approach is that it will fail on a perfectly symmetrical object, but these objects are rare in biological data. \n\nWe apply a series of filters and pick the point based on the first filter it passes:\n\n1. The voxel closest to the centroid of the current label.\n2. The voxel closest to the centroid of the image plane.\n3. Closest to a corner of the plane.\n4. Closest to an edge of the plane.\n5. The previously found maxima.\n\nIt is important that filter #1 be based on the shape of the label so that kinks are minimimized for convex hulls. For example, originally we used only filters two thru five, but this caused skeletons for neurites located away from the center of a chunk to suddenly jink towards the center of the chunk at chunk boundaries.\n\n### Rolling Invalidation Cube\n\nThe original TEASAR paper calls for a \"rolling invalidation ball\" that erases foreground voxels in step 6(iii). A naive implementation of this ball is very expensive as each voxel in the path requires its own ball, and many of these voxels overlap. In some cases, it is possible that the whole volume will need to be pointlessly reevaluated for every voxel along the path from root to target. While it's possible to special case the worst case, in the more common general case, a large amount of duplicate effort is expended.\n\nTherefore, we applied an algorithm using topological cues to perform the invalidation operation in linear time. For simplicity of implmentation, we substituted a cube shape instead of a sphere. The function name `roll_invalidation_cube` is intended to evoke this awkwardness, though it hasn't appeared to have been important. \n\nThe two-pass algorithm is as follows. Given a binary image *I*, a skeleton *S*, and a set of vertices *V*:\n\n1. Let *Bv* be the set of bounding boxes that inscribe the spheres indicated by the TEASAR paper.\n2. Allocate a 3D signed integer array, *T*, the size and dimension of *I* representing the topology. *T* is initially set to all zeros.\n3. For each *Bv*:\n 1. Set T(p) += 1 for all points *p* on *Bv*'s left boundary along the x-axis.\n 2. Set T(p) -= 1 for all points *p* on *Bv*'s right boundary along the x-axis.\n4. Compute the bounding box *Bglobal* that inscribes the union of all *Bv*.\n5. A point *p* travels along the x-axis for each row of *Bglobal* starting on the YZ plane. \n 1. Set integer *coloring* = 0\n 2. At each index, *coloring* += *T*(p)\n 3. If *coloring* > 0 or *T*(p) is non-zero (we're on the leaving edge), we are inside an invalidation cube and start converting foreground voxels into background voxels.\n\n## Related Projects\n\nSeveral classic algorithms had to be specially tuned to make this module possible. \n\n1. [edt](https://github.com/seung-lab/euclidean-distance-transform-3d): A single pass, multi-label anisotropy supporting euclidean distance transform implementation. \n2. [dijkstra3d](https://github.com/seung-lab/dijkstra3d): Dijkstra's shortest-path algorithm defined on 26-connected 3D images. This avoids the time cost of edge generation and wasted memory of a graph representation.\n3. [connected-components-3d](https://github.com/seung-lab/connected-components-3d): A connected components implementation defined on 26-connected 3D images with multiple labels.\n4. [fastremap](https://github.com/seung-lab/fastremap): Allows high speed renumbering of labels from 1 in a 3D array in order to reduce memory consumption caused by unnecessarily large 32 and 64-bit labels.\n\nThis module was originally designed to be used with CloudVolume and Igneous. \n\n1. [CloudVolume](https://github.com/seung-lab/cloud-volume): Serverless client for reading and writing petascale chunked images of neural tissue, meshes, and skeletons.\n2. [Igneous](https://github.com/seung-lab/igneous/tree/master/igneous): Distributed computation for visualizing connectomics datasets. \n\nSome of the TEASAR modifications used in this package were first demonstrated by Alex Bae.\n\n1. [skeletonization](https://github.com/seung-lab/skeletonization): Python implementation of modified TEASAR for sparse labels.\n\n## Credits \n\nAlex Bae developed the precursor skeletonization package and several modifications to TEASAR that we use in this package. Alex also developed the postprocessing approach used for stitching skeletons using 50% overlap. Will Silversmith adapted these techniques for mass production, refined several basic algorithms for handling thousands of labels at once, and rewrote them into the Kimimaro package. Will added trickle DAF, zero weighted previously explored paths, and fixing borders to the algorithm. Forrest Collman added parameter flexibility and helped tune DAF computation performance. Sven Dorkenwald and Forrest both provided helpful discussions and feedback.\n\n## References \n\n1. M. Sato, I. Bitter, M.A. Bender, A.E. Kaufman, and M. Nakajima. \"TEASAR: Tree-structure Extraction Algorithm for Accurate and Robust Skeletons\". Proc. 8th Pacific Conf. on Computer Graphics and Applications. Oct. 2000. doi: 10.1109/PCCGA.2000.883951 ([link](https://ieeexplore.ieee.org/abstract/document/883951/))\n2. I. Bitter, A.E. Kaufman, and M. Sato. \"Penalized-distance volumetric skeleton algorithm\". IEEE Transactions on Visualization and Computer Graphics Vol. 7, Iss. 3, Jul-Sep 2001. doi: 10.1109/2945.942688 ([link](https://ieeexplore.ieee.org/abstract/document/942688/))\n3. T. Zhao, S. Plaza. \"Automatic Neuron Type Identification by Neurite Localization in the Drosophila Medulla\". Sept. 2014. arXiv:1409.1892 \\[q-bio.NC\\] ([link](https://arxiv.org/abs/1409.1892))\n4. A. Tagliasacchi, T. Delame, M. Spagnuolo, N. Amenta, A. Telea. \"3D Skeletons: A State-of-the-Art Report\". May 2016. Computer Graphics Forum. Vol. 35, Iss. 2. doi: 10.1111/cgf.12865 ([link](https://onlinelibrary.wiley.com/doi/full/10.1111/cgf.12865))\n5. P. Li, L. Lindsey, M. Januszewski, Z. Zheng, A. Bates, I. Taisz, M. Tyka, M. Nichols, F. Li, E. Perlman, J. Maitin-Shepard, T. Blakely, L. Leavitt, G. Jefferis, D. Bock, V. Jain. \"Automated Reconstruction of a Serial-Section EM \u200b Drosophila \u200b Brain with Flood-Filling Networks and Local Realignment\". April 2019. bioRXiv. doi: 10.1101/605634 ([link](https://www.biorxiv.org/content/10.1101/605634v1))\n6. M.M. McKerns, L. Strand, T. Sullivan, A. Fang, M.A.G. Aivazis, \"Building a framework for predictive science\", Proceedings of the 10th Python in Science Conference, 2011; http://arxiv.org/pdf/1202.1056\n7. Michael McKerns and Michael Aivazis, \"pathos: a framework for heterogeneous computing\", 2010- ; http://trac.mystic.cacr.caltech.edu/project/pathos",
"description_content_type": "text/markdown",
"docs_url": null,
"download_url": "",
"downloads": {
"last_day": -1,
"last_month": -1,
"last_week": -1
},
"home_page": "https://github.com/seung-lab/kimimaro/",
"keywords": "",
"license": "GPLv3+",
"maintainer": "",
"maintainer_email": "",
"name": "kimimaro",
"package_url": "https://pypi.org/project/kimimaro/",
"platform": "",
"project_url": "https://pypi.org/project/kimimaro/",
"project_urls": {
"Homepage": "https://github.com/seung-lab/kimimaro/"
},
"release_url": "https://pypi.org/project/kimimaro/1.0.1/",
"requires_dist": null,
"requires_python": "",
"summary": "Skeletonize densely labeled image volumes.",
"version": "1.0.1"
},
"last_serial": 5986253,
"releases": {
"0.1.0": [
{
"comment_text": "",
"digests": {
"md5": "0e12328535ab5aad8bbab6e25eee6c6a",
"sha256": "b9daf71907121654315619ddd5c6099722607dc5afcc84b295e6730ac45f0eaa"
},
"downloads": -1,
"filename": "kimimaro-0.1.0-cp36-cp36m-macosx_10_13_x86_64.whl",
"has_sig": false,
"md5_digest": "0e12328535ab5aad8bbab6e25eee6c6a",
"packagetype": "bdist_wheel",
"python_version": "cp36",
"requires_python": null,
"size": 133927,
"upload_time": "2019-01-28T04:20:04",
"url": "https://files.pythonhosted.org/packages/fb/d1/7383f939c73be4b980111fc0191951c35f6a50dc74fc2196c384209eaa6f/kimimaro-0.1.0-cp36-cp36m-macosx_10_13_x86_64.whl"
},
{
"comment_text": "",
"digests": {
"md5": "35a51fa4876e17d3bee5a8fa1b5caa57",
"sha256": "aa5614564724af29f0a6964e93fac956623887b8bdc075f617ce061bb3e2b0fb"
},
"downloads": -1,
"filename": "kimimaro-0.1.0.tar.gz",
"has_sig": false,
"md5_digest": "35a51fa4876e17d3bee5a8fa1b5caa57",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 631801,
"upload_time": "2019-01-28T04:20:07",
"url": "https://files.pythonhosted.org/packages/42/32/fec8b362789028ec1d738c16eb6ccb5635533040a6e0307f108f0dce8249/kimimaro-0.1.0.tar.gz"
}
],
"0.2.0": [
{
"comment_text": "",
"digests": {
"md5": "55b65115b0e67920698b4b65b12ba890",
"sha256": "49d893bd92e8a52f47b6a14eb246ef612da2290e39246e7718f036697ce933d4"
},
"downloads": -1,
"filename": "kimimaro-0.2.0.tar.gz",
"has_sig": false,
"md5_digest": "55b65115b0e67920698b4b65b12ba890",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 633569,
"upload_time": "2019-01-30T03:59:27",
"url": "https://files.pythonhosted.org/packages/df/8d/130759db41829cfaf963c99afa80e710a1e59089d661485ea82e7f15abdb/kimimaro-0.2.0.tar.gz"
}
],
"0.2.1": [
{
"comment_text": "",
"digests": {
"md5": "f18344dfa93477af981c2f7577c9ec16",
"sha256": "2c0b05b2ca2e9227aac21b7a17ce81ba45eeea28e3d251890c347a904a0098a6"
},
"downloads": -1,
"filename": "kimimaro-0.2.1.tar.gz",
"has_sig": false,
"md5_digest": "f18344dfa93477af981c2f7577c9ec16",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 633592,
"upload_time": "2019-02-01T23:32:55",
"url": "https://files.pythonhosted.org/packages/0a/8f/18a802c35f4a6769437eb013f9dc2bbeaf300800579f22447093cbfdf0c6/kimimaro-0.2.1.tar.gz"
}
],
"0.2.2": [
{
"comment_text": "",
"digests": {
"md5": "3a54a82993f9021273b0a99c0ad7a108",
"sha256": "e3f89fccc048625144deb776d32f5324f26a88eecab1b83d1a06f807539c2a1a"
},
"downloads": -1,
"filename": "kimimaro-0.2.2.tar.gz",
"has_sig": false,
"md5_digest": "3a54a82993f9021273b0a99c0ad7a108",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 875438,
"upload_time": "2019-03-19T22:47:29",
"url": "https://files.pythonhosted.org/packages/bf/8b/26c1c49079315f1fe126780df63e2bdad2dfde6192e1e9220d18da333dee/kimimaro-0.2.2.tar.gz"
}
],
"0.3.0": [
{
"comment_text": "",
"digests": {
"md5": "4ba78b88adcbabb4d14b6392447ed99a",
"sha256": "4d84b1f1b93c1a4ebfea8a590e7d19f83c207632fdb51a0d17200f36d26864a5"
},
"downloads": -1,
"filename": "kimimaro-0.3.0.tar.gz",
"has_sig": false,
"md5_digest": "4ba78b88adcbabb4d14b6392447ed99a",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 903559,
"upload_time": "2019-04-25T01:40:17",
"url": "https://files.pythonhosted.org/packages/7f/ba/ab3db0d96b241c038d8e11b864a87a5b70137e5a86938beb3bda9c43f849/kimimaro-0.3.0.tar.gz"
}
],
"0.3.1": [
{
"comment_text": "",
"digests": {
"md5": "21c5fceab43f3dbb575a231b7aa730a3",
"sha256": "fc4b479c51dd56199ffe2b7dbde6252e7ba9fc524ae9536af20dae46c56034fb"
},
"downloads": -1,
"filename": "kimimaro-0.3.1.tar.gz",
"has_sig": false,
"md5_digest": "21c5fceab43f3dbb575a231b7aa730a3",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 911355,
"upload_time": "2019-04-25T04:11:34",
"url": "https://files.pythonhosted.org/packages/5c/13/0889ec9c1387ad233200ddff719a44fd4f11a71917a7a0259abc3a86210b/kimimaro-0.3.1.tar.gz"
}
],
"0.4.0": [
{
"comment_text": "",
"digests": {
"md5": "646ad3dc2456390b083c6237ab76258a",
"sha256": "76b465327fa8fdc8ba84b9ac389ea99ed5966b52b3725630c84c2d1ddd93b038"
},
"downloads": -1,
"filename": "kimimaro-0.4.0.tar.gz",
"has_sig": false,
"md5_digest": "646ad3dc2456390b083c6237ab76258a",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 912644,
"upload_time": "2019-04-26T20:03:47",
"url": "https://files.pythonhosted.org/packages/8e/d7/6854957ce99642e0cb99760042da6ed0afc144d219d44ddc5123bd791cd4/kimimaro-0.4.0.tar.gz"
}
],
"0.4.1": [
{
"comment_text": "",
"digests": {
"md5": "8b2029d0ae5bf2c377cb3aaf6d0b9ea1",
"sha256": "c78015e6c018d66f172f6ac509a08bccaa31718b61f9787d9492fd50a8be4eb7"
},
"downloads": -1,
"filename": "kimimaro-0.4.1.tar.gz",
"has_sig": false,
"md5_digest": "8b2029d0ae5bf2c377cb3aaf6d0b9ea1",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 915273,
"upload_time": "2019-05-08T00:50:05",
"url": "https://files.pythonhosted.org/packages/57/b5/0ac2023a096fec784c23a999bc93ec1339b357a4a20375bd973da36f0505/kimimaro-0.4.1.tar.gz"
}
],
"0.4.2": [
{
"comment_text": "",
"digests": {
"md5": "381545ca4d7ff75f998d09e93e44a092",
"sha256": "b6605c2e50a5f05e6f8ab3cb4c2be135438c5195ff877087a8789e36d7d8c9e7"
},
"downloads": -1,
"filename": "kimimaro-0.4.2.tar.gz",
"has_sig": false,
"md5_digest": "381545ca4d7ff75f998d09e93e44a092",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 915482,
"upload_time": "2019-05-30T03:53:31",
"url": "https://files.pythonhosted.org/packages/79/d0/8cb85e2ae99e5651578db83850850a18e0102287593bcb5a2f545adfdd3e/kimimaro-0.4.2.tar.gz"
}
],
"0.5.0": [
{
"comment_text": "",
"digests": {
"md5": "5fc6520f218a2ab640c10a02b512f7d6",
"sha256": "f3e51f66b225bc17fa3bea2fcb056bb01aebaf7f8030ddebca602278386981e5"
},
"downloads": -1,
"filename": "kimimaro-0.5.0.tar.gz",
"has_sig": false,
"md5_digest": "5fc6520f218a2ab640c10a02b512f7d6",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1096038,
"upload_time": "2019-06-13T18:14:44",
"url": "https://files.pythonhosted.org/packages/f1/90/14ffa60f1a6f6b18c0b1af91de41319b50d8f879b85843c00cdd3258ed19/kimimaro-0.5.0.tar.gz"
}
],
"0.5.1": [
{
"comment_text": "",
"digests": {
"md5": "4740755298cbca2923f754c83c02bb25",
"sha256": "85c3135b7bce1805adfc9702e8dd3fc800d1a11e80e303c622da2a2c67d62821"
},
"downloads": -1,
"filename": "kimimaro-0.5.1.tar.gz",
"has_sig": false,
"md5_digest": "4740755298cbca2923f754c83c02bb25",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1108508,
"upload_time": "2019-06-16T01:11:08",
"url": "https://files.pythonhosted.org/packages/c0/8e/c731a9eecd2d5fc3e198ba5a677bc1c3ca3057a1ff61fd854a2a34743ae6/kimimaro-0.5.1.tar.gz"
}
],
"0.5.2": [
{
"comment_text": "",
"digests": {
"md5": "340da76a03fc902cb6df0c458bbd26fb",
"sha256": "16fc21327aaaf20102ed686b8028641a42a8d8578d0b356dc841810aaf94073c"
},
"downloads": -1,
"filename": "kimimaro-0.5.2.tar.gz",
"has_sig": false,
"md5_digest": "340da76a03fc902cb6df0c458bbd26fb",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1108661,
"upload_time": "2019-06-21T18:18:40",
"url": "https://files.pythonhosted.org/packages/71/e4/9b6d71103c1bf0280de3ab73bff5dc1f0b8f4f8ce1dc0186fdd2c36904e4/kimimaro-0.5.2.tar.gz"
}
],
"0.5.3": [
{
"comment_text": "",
"digests": {
"md5": "736d3fcc9ba5ac99ce6e3da6c44b563e",
"sha256": "3662a8b75bfe88db545293a7126670c41b3c26ca2cad28b8810e974e452313fa"
},
"downloads": -1,
"filename": "kimimaro-0.5.3.tar.gz",
"has_sig": false,
"md5_digest": "736d3fcc9ba5ac99ce6e3da6c44b563e",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1108706,
"upload_time": "2019-07-25T22:59:49",
"url": "https://files.pythonhosted.org/packages/18/7e/5ec11e137bc585bff692d56ae51034b902565ae1cdedde95efd5cbd60390/kimimaro-0.5.3.tar.gz"
}
],
"0.5.4": [
{
"comment_text": "",
"digests": {
"md5": "3db2b82b19320356321d0e039e3da092",
"sha256": "da34d467246862fff54b1a2a44875e2baf58ff97242b8aefa8525bc82767eaf6"
},
"downloads": -1,
"filename": "kimimaro-0.5.4.tar.gz",
"has_sig": false,
"md5_digest": "3db2b82b19320356321d0e039e3da092",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1268222,
"upload_time": "2019-08-02T21:21:44",
"url": "https://files.pythonhosted.org/packages/d1/d1/292e38f6e32c7942258c9b1a22c60ac05a101a6503947efaf3c59d85429a/kimimaro-0.5.4.tar.gz"
}
],
"0.6.0": [
{
"comment_text": "",
"digests": {
"md5": "28f9ced56285a00d7c57dfb413642b77",
"sha256": "b74864bc671435b1f3b42bc517de1d01226e90deea4a6be7d6565c8aa4cc4445"
},
"downloads": -1,
"filename": "kimimaro-0.6.0.tar.gz",
"has_sig": false,
"md5_digest": "28f9ced56285a00d7c57dfb413642b77",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1274282,
"upload_time": "2019-08-20T21:14:00",
"url": "https://files.pythonhosted.org/packages/e7/6b/75f17720e2cb7483f441c79a9b117a4c04e808e2f3f2e05e0dc43f01dba9/kimimaro-0.6.0.tar.gz"
}
],
"0.7.0": [
{
"comment_text": "",
"digests": {
"md5": "41e51f778a7ea1b669770639b11b4c63",
"sha256": "259385c68847e045c50f2188d970bb9161e323c12cd4e1d497b7134027d22c56"
},
"downloads": -1,
"filename": "kimimaro-0.7.0.tar.gz",
"has_sig": false,
"md5_digest": "41e51f778a7ea1b669770639b11b4c63",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1276225,
"upload_time": "2019-09-10T05:31:14",
"url": "https://files.pythonhosted.org/packages/16/6a/b8368cd2356b75906d13dbd5e7b11061089da9010844eda178e5d2e120e4/kimimaro-0.7.0.tar.gz"
}
],
"1.0.0": [
{
"comment_text": "",
"digests": {
"md5": "ab6d2d4308a70cff9aeaffc30fe1a8d2",
"sha256": "bc2e44af33980143fede645fc07005a82cf4900ab85f9eed212f0f32702608ff"
},
"downloads": -1,
"filename": "kimimaro-1.0.0.tar.gz",
"has_sig": false,
"md5_digest": "ab6d2d4308a70cff9aeaffc30fe1a8d2",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1277621,
"upload_time": "2019-10-14T22:49:47",
"url": "https://files.pythonhosted.org/packages/ab/e1/da952cbae989d5844f78c7a8c6332e2040fb6a96e0e2002c3fbe238a0d35/kimimaro-1.0.0.tar.gz"
}
],
"1.0.1": [
{
"comment_text": "",
"digests": {
"md5": "d7db9103885543084926ffba525cb925",
"sha256": "c6e45620320c54650c505a203470d7e380c681980a28dbab85c594544868b016"
},
"downloads": -1,
"filename": "kimimaro-1.0.1.tar.gz",
"has_sig": false,
"md5_digest": "d7db9103885543084926ffba525cb925",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1277725,
"upload_time": "2019-10-16T20:28:50",
"url": "https://files.pythonhosted.org/packages/f4/d6/d193a0848e940c0cae213b789ee7e81b034e0ac1a6ba90158025bfe911f5/kimimaro-1.0.1.tar.gz"
}
]
},
"urls": [
{
"comment_text": "",
"digests": {
"md5": "d7db9103885543084926ffba525cb925",
"sha256": "c6e45620320c54650c505a203470d7e380c681980a28dbab85c594544868b016"
},
"downloads": -1,
"filename": "kimimaro-1.0.1.tar.gz",
"has_sig": false,
"md5_digest": "d7db9103885543084926ffba525cb925",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 1277725,
"upload_time": "2019-10-16T20:28:50",
"url": "https://files.pythonhosted.org/packages/f4/d6/d193a0848e940c0cae213b789ee7e81b034e0ac1a6ba90158025bfe911f5/kimimaro-1.0.1.tar.gz"
}
]
}
![\"A](\"https://raw.githubusercontent.com/seung-lab/kimimaro/master/mass_skeletonization.png\")
![\"A](\"https://raw.githubusercontent.com/seung-lab/kimimaro/master/kimimaro_512x512x512_benchmark.png\")