Many labs have developed separate solutions for storing, visualizing, and annotating datasets of ever-increasing size as contemporary commercial solutions were not adequately scalable, were missing features, or both. These tools operate on many different principles, and most of them have a method for importing images. However, there is limited support for meshing and skeletonization of dense segmentation in bulk as we will describe below.

To give a brief sketch of the landscape, these systems can be broadly characterized by the maximum data size they support (in-memory, disk, network filesystem / cloud storage), method of neural representation (skeletons, per-voxel image segmentation, network representation), their methods of visualization [single resolution image, image pyramid (Pietzsch et al., 2015; Sofroniew et al., 2022), precomputed or dynamic surface meshes, skeletons, volume rendering (Peng et al., 2010; Maitin-Shepard et al., 2021), or other], method of data storage [single file, chunks, random-access consolidated files (“shards”), versioning, and serverlessness], method of proofreading [skeleton tracing (Saalfeld et al., 2009; Boergens et al., 2017) and segment merging and splitting (Kim et al., 2014; Katz and Plaza, 2019; Dorkenwald et al., 2021), with variations within each category], and openness of access to reconstructions and proofreading prior to publication (public, semi-public, or internal only), and the different overlapping communities of proofreaders, viewers, and software developers centered around each tool. Proofreading systems for editing the annotations and segmentation may have additional storage requirements such as a graph of connected segments (Anderson et al., 2011; Ai-Awami et al., 2016; Dorkenwald et al., 2021). An excellent overview of the visualization landscape can be found in a recent survey by Beyer et al. (2022) and in the VAST paper by Berger et al. (2018).

To place Igneous in relation to these categories, it produces Neuroglancer Precomputed volumes that are stored in cloud or network storage. Neurons are represented by per-voxel annotations, precomputed multi-resolution surface meshes, and precomputed skeletons. The image data are stored as versionless chunks or shards (see Consolidating Files). Igneous is free software that can be used by anyone. See the Supplementary material for additional information.

Random read access to large datasets is often achieved by chunking large images into many smaller files or writing out each mesh, skeleton, or other derivatives of a label as one or more files. This results in problems for the filesystem when the number of files becomes large. Even with cloud storage, as of this writing, file creation costs between $5²⁶ and $6.50²⁷ per a million files on several cloud providers on active storage (usually the default tier). When hundreds of millions or billions of files are created, the initial upload costs can be more than the cost of the computation.

While it's difficult to impute the financial cost directly to the stress put on the storage system, it is notable that object storage systems often use an erasure coding (Balaji et al., 2018) or replication scheme that creates multiple copies or parity fragments to be able to reconstruct each file in the event of bit flips or damage to machines and hard disks. Frequently, the metadata used to locate distributed copies or file fragments are themselves replicated several times. Thus, there is a constant storage multiplier attached to each file uploaded. Sometimes, even if a large dataset can be created, deleting it can be difficult on some systems, incurring delays or additional costs.

As a mitigation, Neuroglancer adopted an approach that condenses many individual files into random-read (but not random-write) files called “shards”²⁸ (see Figure 4). Shards reduce the number of files in a dataset by several orders of magnitude. This kind of approach is becoming more popular and implementations are being discussed in prominent projects like Zarr.²⁹ Techniques similar to sharding can be implemented by a filesystem to cope with large numbers of small files, so shards could be viewed as a client-side emulation of a filesystem feature.

Igneous provides methods for creating shard files for each data type it supports (images, meshes, and skeletons) and CloudVolume is capable of reading them. As producing shards is more complicated and removes random access writes, it is usually hidden behind a flag (e.g., --sharded) or a separate command (e.g., merge-sharded).

Igneous supports all encoding methods that are currently supported in Neuroglancer Precomputed.

Large volumetric images do not fit in RAM and are not practical to download for visualization. It is also desirable to perform some types of processing on lower resolution images. Thus, it is necessary to create an image pyramid of lower resolution images that can be downloaded selectively (e.g., depending on the viewer's current zoom level). Typically, each successive resolution level (or “mip”⁴⁶) is downsampled using 2 × 2 (2D) or 2 × 2 × 2 (3D) voxel pooling. For Neuroglancer, mip 0 is the highest resolution layer, with higher integers indicating successively lower resolution layers.

In calculating these mip levels, we use average pooling for natural (electron microscopy) images and recursive mode pooling for labeled (segmentation) images. Uncompressed, a pyramid of 2 × 2 downsamples requires 33% additional storage over the base image. A pyramid of 2 × 2 × 2 downsamples will require 14% more. In both cases, mip 1 is the dominant contributor to the additional storage (25% for 2 × 2, 12.5% for 2 × 2 × 2).

Multiple mip levels can be generated from a single task if the task is large enough. This is advantageous because (a) it takes fewer runs to generate all mip levels (b) fewer tasks need to be generated and managed (c) the average pooling algorithm can avoid integer truncation. Nonetheless, it's often not possible to generate all mip levels in a single task as each additional level exponentially multiplies the required memory by 4–8 times. Therefore, by default, five levels are generated and another set can be generated on top of them in a process we term “superdownsampling.”

Integer truncation becomes significant in very large volumes that may have 10+ mip levels. If it is not accounted for, the lowest resolution levels are noticeably darker than the base level as up to 0.75 luminance units can be lost per voxel at each level in 2 × 2 pooling or 0.875 units in 2 × 2 × 2 pooling. At 10 mips, this amounts to a potential loss of 7.5 units in 2 × 2 and 8.75 for 2 × 2 × 2. Our average pooling implementation in the tinybrain⁴⁷ library mitigates this issue by accumulating sums for each mip level before dividing.

2 × 2 × 2 downsampling can lead to ghosting around the edges of a slice as extensions of one slice regions can be averaged with several zeros from a gap in an adjacent slice. We offer a sparse downsampling mode that skips counting zero-valued voxels in the average.

For sharded volumes, due to the large memory requirement for holding a single shard in memory, Igneous can only generate a single mip level at a time. Generating additional mip levels would require holding exponentially larger numbers of shards in memory for each mip level unless the shards shrink at each mip level, reducing their utility. This requirement means that producing a sharded mip incurs integer truncation errors. However, it is possible to produce unsharded mips from sharded mips and later condense them.

While visualizing cross sections of 3D segmentation are informative, it is difficult to understand the overall shape of a neuron outside of a series of local contexts without a visual representation of the whole object. Meshing, creating a 3D surface representation consisting of vertices and faces from a segmentation, provides a light-weight method of visualizing neurons that may span an entire dataset as only the surface of the object of interest needs to be represented.

Neuroglancer offers three different formats for representing meshes: single-resolution unsharded, multi-resolution unsharded, and multi-resolution sharded (single-resolution sharded is not supported). The single resolution format allows for multiple mesh fragments for each segmentation label to be written as separate files and then pointed to by a manifest JSON file with a well-known filename (“SEGID:0”). The multi-resolution formats are different in that each mesh data file must contain the whole mesh and potentially additional lower resolution meshes addressable in an octree format. In the unsharded version, each label has a data file and an index file that gives instructions for reading the data file's octree. In the sharded format, a mesh's index file and data file are concatenated and grouped with many other labels's meshes. The benefit of multi-resolution files is that Neuroglancer can display lower resolution meshes when the viewer is zoomed out, which can allow for a higher performance display of large or numerous meshes.

We produce meshes by applying a specialized variant of the marching cubes algorithm (Lorensen and Cline, 1987) to segmentation images and then simplify the results. Two passes are needed. As it is impossible to produce a mesh spanning large datasets within the memory of a single machine, the image must be divided into smaller tasks and then merged. The first pass produces mesh fragments derived from small spatial areas and writes them to storage. The second pass aggregates the mesh fragments to create the final output (and in the case of multi-resolution meshes, also generates a hierarchy of simplified meshes).

Typically we begin with downsampling the segmentation image to use a near-isotropic mip level for meshing. While any mip level can be used, downsampling vastly reduces the amount of data to be processed, often by 64 × , while retaining a reasonable amount of detail. The selected mip level's image is then divided into a regular grid of tasks that each overlap on their positive axis face by one voxel.

We wrote the free software zmesh⁴⁸ library to produce our meshes. zmesh is based on zi_lib⁴⁹, which was first used in the Omni (Shearer, 2009) neural reconstruction proofreading software and then in Eyewire (Kim et al., 2014), an online crowdsourced proofreading platform. It implements a specialized version of the marching cubes algorithm (Lorensen and Cline, 1987) which efficiently handles multi-labeled images. This allows for all meshes to be generated in a single pass instead of thousands. The resultant mesh fragments are then simplified using the methods of Garland and Heckbert (1997) and Hoppe (1999) modified to retain topological integrity by filtering out unsound simplifications.

The one voxel overlap causes marching cubes to output mesh fragments that can be trivially stitched. Each fragment is then simplified before being serialized and compressed with gzip. Single resolution meshes are serialized in the legacy Precomputed format. Meshes written to the multi-resolution format are Draco compressed with integer position attributes using DracoPy.⁵⁰ To generate a resolution hierarchy, we repeatedly apply pyfqmr⁵¹ to the base mesh with increasingly aggressive settings.

Skeletons are stick-figure centerline representations of neurons. They are often represented as trees or graphs with vertices that are localized in space. They have many uses in analyzing the anatomy and function of neurons. For example, they can be used to compute basic properties such as cable length and wire diameter (when the skeleton is a medial axis transform). They can also be used to define functional sections of neurite during analysis such as dendrite and axon or to model electrical compartments. Skeletons can be used to guide cameras in 3D renderings and in proofreading neural circuits (as in CATMAID and WebKnossos).

However, mass-producing large skeletons from dense segmentation images is difficult due to both performance and memory requirements. Most image skeletonization algorithms require holding the entire image in memory, which is impossible at near petavoxel scale. The most popular skeletonization algorithm in the connectomics field is the TEASAR algorithm (Sato et al., 2000; Bitter et al., 2001), due to its flexibility in skeletonization different shapes with parameterization (Zhao and Plaza, 2014). However, the previously available implementations of TEASAR are very slow and require high memory usage. Furthermore, like meshes, mass production of skeletons can produce hundreds of millions of files, which results in difficulties with the filesystem or high cost on cloud storage. These problems previously made the mass production of skeletons impractical. However, we have been able to overcome these limitations and can now mass produce skeletons directly from segmentation images and without the need for a meshing step.

We developed Kimimaro (Silversmith et al., 2021a), a fast Python, Cython, and C++ based TEASAR-like algorithm that can process densely labeled segmentation images at hundreds of kilovoxels per a second with memory usage low enough to use on consumer hardware with a 512³ voxel cutout. However, even with improvements in the core algorithm, it was still not feasible to skeletonize the entire image in one shot due to memory constraints. Thus, we chunk the image into a regular grid of 512³ voxel tasks and stitch the results together in a second pass, similar to meshing. For this approach to work, several conditions must be satisfied (a) the stitching process must be made reliable (in the sense that the right connection between tasks is always made) and efficient (b) the tasks must be large enough to encompass significant morphological features or else skeletonization may go awry from lack of context (c) stitched skeletons should be free of loops.

Igneous ensures reliable connections between skeleton chunks by using single voxel overlap between tasks and ensuring that skeletons generated on each side of the border will meet at the same voxel. This property allows them to be trivially stitched before post-processing. This is accomplished by changing the TEASAR algorithm to first target borders before proceeding normally. A border target is defined for each 2D connected component extracted from the boundary cross-section of each shape. Each target is defined as the voxel containing the peak euclidean distance transform of each connected component (see Figure 5). This metric ensures that, unlike a centroid, the target voxel resides within the component.

To break ties between peak voxels, we use topological features to ensure that the criteria are coordinate-frame independent as the connecting coordinate frames are mirrored on each side of the task (top & bottom, left & right, back & front). The tie-breaking features in order of precedence:

These criteria work well for most shapes but a small number of shapes such as an annulus, X, square shape, or specially constructed irregular shapes centered at the image center will generate up to eight candidate points. We haven't implemented a final tie-breaker, but it could be resolved by introducing an asymmetric criterion (such as selecting the top-left peak). This can be made to work if it alternates between top-left and top-right depending on the orientation (e.g., up vs. down) of the face, otherwise, it will introduce a disconnection.

If a shape contacts the edge or corner of the bounding box, it will generate two and three targets respectively. The adjacent tasks will also draw skeletons to each of these target points, resulting in small loops. These loops are later resolved during post-processing during which loops are removed, small extensions are pruned, and components closer together than the nearest boundary are joined.

By using a fast and memory-efficient library, a chunked skeletonization strategy, and supporting sharded skeleton production (see Condensing Files), Igneous makes mass production of quality skeletons practical. Ensuring that Kimimaro has appropriate visual context is a somewhat more difficult problem that doesn't yet have a perfect solution as the arbitrary division of the image into a regular grid can split objects. A reasonable way to manage this problem is to ensure that the amount of visual context in the image is greater than the size of the largest objects of interest. Reducing image resolution and increasing the size of each task can both increase visual context.

Similar to meshing, we often skeletonize volumes at a near-isotropic resolution. While our highest resolution segmentation is usually 4 × 4 × 40 or 8 × 8 × 40 nm³, ordinarily skeletonization is processed at 32 × 32 × 40 nm³. At lower resolutions, such as 64 × 64 × 40 nm³, we find that skeletons become noticeably de-centered from neurites. At even lower resolutions, thin processes may become disconnected.

By default, skeletonization tasks are 512³ voxels, which at 32 × 32 × 40 nm³ equates to 16.4 × 16.4 × 20.5 μm³ of physical context. The interconnection scheme described above works well for wire-like objects, but more bulbous objects such as somata, require full context to produce a reasonable skeleton. It is difficult to guarantee that large objects will have full context when crudely dividing the image, so future work will be needed to refine this aspect.

Contrast correction via histogram equalization on a per-Z slice basis is supported for 8 and 16-bit images with optional right and left tail clipping. This operation is two pass as statistical information about the histogram of each slice must be known before the adjustment can be made. The first pass collects sample data from patches assigned in a regular grid at a configurable fraction of the image area (default 1%). It then writes a JSON file for each slice containing the sampled histogram. A second pass of contrast correction tasks then performs the histogram equalization on a regular grid of chunks to avoid downloading an entire slice (which may be dozens of gigabytes). Downsampling is integrated into the second pass to enable rapid visualization and save work.

Managing large datasets is a problem in its own right. Simply enumerating a billion files can become a challenge when only thousands or tens of thousands of filenames can be listed per second. Nonetheless, it is frequently advantageous to move datasets to share them with collaborators, to move the data closer to the site of computation, to use a more economical storage provider, or to re-encode them with a different compression algorithm.

Igneous provides convenient commands to transfer, re-encode, re-chunk, delete, and condense large data for images, meshes, and skeletons. Transfers and deletions are often more efficient when working with sharded volumes as fewer files need to be manipulated.

Igneous can be installed on any system meeting the following requirements (a) runs a supported cPython version (currently 3.7+) (b) runs a recent Linux, MacOS, or Windows operating system OR can run an Ubuntu Linux based Docker container (c) all workers can access a common queue via either FileQueue which requires a filesystem with consistent advisory file locks OR via internet access to AWS SQS. Installation with Python pip is simple: pip install igneous-pipeline.

Igneous can be scripted in Python or run from the command line. A typical workflow is to first select a dataset and operation, and then enqueue a set of tasks in either FileQueue or Amazon SQS from a local workstation. Then, using the execute command, which may be used on the same workstation or a cluster, the queue is selected and executed against with one worker process per an available core. Execution continues until the queue is empty. Termination can be set to automatic in the case of FileQueue, but SQS only returns the approximate number of tasks enqueued and so requires either manual monitoring or repeated polling to verify the current job is finished. Examples of Igneous CLI commands can be seen in Listing 1 and an enumeration of the available commands is available in Table 1.

Parallel execution can be accomplished by any runner that can run the command igneous execute$QUEUE where $QUEUE is a path (such as ./queue or sqs://my-queue). On a local machine, parallel execution can be triggered by the -p flag which indicates the number of processes to start. For large distributed jobs, we have used Docker/Kubernetes and SLURM as runners, but many other platforms would be suitable. Queue progress can be monitored and managed via the co-installed ptq (short for “Python Task Queue”) command line utility.

As images generally comprise the bulk of storage for a given connectomics dataset, we attempted to characterize the efficacy of different compression technologies supported by Neuroglancer. We downloaded three pairs of padded test image and training segmentation datasets. We converted the HDF5 files into uncompressed Neuroglancer Precomputed volumes using CloudVolume and then re-encoded the raw volume using the igneous xfer command. For each re-encoding, we measured the computation time taken and the resultant disk space used. For decoding, we timed reading each resultant volume with CloudVolume 8.7.0 five times and reported the average.

All measurements were taken on a 2021 M1 Macbook Pro with an SSD using one process. The measurements for CREMI volumes A+, B+, and C+ were averaged into a point estimate of computation time in megavoxels per second (MVx/s) and disk space required for each encoding. All three CREMI image volumes were 8-bit 3072 × 3072 × 200 voxels (1.9 gigavoxels each, 5.7 gigavoxels total). All three CREMI segmentation volumes were all rendered as 64-bit 1250 × 1250 × 125 voxels (195 MVx each, 586 MVx total).

For microscopy images (see Figure 6), raw (uncompressed), raw+gzip, raw+brotli, jpeg, and png transfer encodings were evaluated. PNG and gzip were evaluated at compression level 9, brotli at level 5, and JPEG at 85% quality. For segmentation images, we evaluated raw (uncompressed), compresso, and compressed_segmentation encodings layered in combination with gzip level 9 and brotli level 5.

The best compression and encoding speed for images was the lossy JPEG encoding (5.98 × , 288.6 MVx/s). The best lossless compression was PNG (1.89 × ) though it was the slowest (13.7 MVx/s). raw+brotli and raw+gzip had similar encoding speeds (25.7 and 25.8 MVx/s), but raw+brotli gave a slightly smaller file size (1.40 × gzip and 1.48 × brotli).

Decoding was uniformly much faster than encoding for microscopy images. The fastest decoding time was naturally held by uncompressed data (753.4 MVx/s). jpeg was second (574.5 MVx/s). Of the lossless compression codecs, raw+gzip (363.0 MVx/s) and raw+brotli (361.8 MVx/s) were similar. png was the slowest (74.5 MVx/s).

For segmentation (see Figure 7) we evaluated raw, compressed_segmentation, and compresso each with no compression, gzip, and brotli second stage compression. The overall best compression was given by compresso+brotli (570 × , 64.7 MVx/s). compresso+gzip gave a similar compression ratio, but was slightly slower at encoding. compressed_segmentation+brotli was slightly faster (69.4 MVx/s) but only yielded a 333 × compression ratio (58%). The fastest overall method was writing uncompressed raw arrays (98.5 MVx/s), though it is not suitable for large datasets. Using compresso+brotli resulted in an 8 × speedup and a 4.8 × compression improvement vs. the naïve raw+gzip approach.

Decoding speeds were faster than encoding speeds in all segmentation trials and showed less variability between codecs. Encoding and decoding speeds were also closer in magnitude than with microscopy image compression. compresso was slowest (96–97 MVx/s). Interestingly, it seemed to yield the same decoding speed regardless of second layer compression type (including none). raw compression was slightly faster and compressed_segmentation was the fastest encoding type (111–115 MVx/s). raw without any compression was fastest at 122.5 MVx/s.

In order to characterize the computation and disk space required to process a dataset using Igneous, we processed a rough automatic segmentation of the well-known mouse S1 dataset (Kasthuri et al., 2015) which is 283.1 gigavoxels (7.4 GB gzipped) and has a resolution of 6 × 6 × 30 nm³ (see Figure 1). We downloaded the dataset from cloud storage to an 8-core 2021 M1 Macbook Pro with an SSD. The segmentation was then converted into a sharded compresso+gzip encoded image (shards do not currently support brotli) with 256 × 256 × 32 voxel chunks (Compresso is more effective on larger XY planes, hence the asymmetry). It was then unsharded downsampled to mip 5 in one step. The downsampled chunks were compresso+brotli compressed. It was then meshed in Draco compressed sharded format with a single resolution at image resolution 24 × 24 × 30 nm³, and skeletonized using parameters intended to capture spines at the same resolution. The results can be seen in Figure 8.

Table 2 shows the time and disk usage required by each stage of processing. In total, 36 core-hours were required, which works out to 7.9 gigavoxels/core-hr. The most computationally demanding tasks by far were generating mesh and skeleton fragments (“forging”) at all spatial grid points. However, in our experience, if parameters are poorly chosen or large mergers are present, skeleton merging can become demanding as well.

Quite sensibly, segmentation images including downsamples comprise the bulk of disk storage (5.4 GB), but thanks to compression, the highest resolution image (2.3 GB) isn't the single largest consumer of disk space. Instead, mesh intermediate files (4.6 GB) are the largest single item. While the intermediate fragments are simplified, they are gzip compressed while the final meshes will be merged and Draco compressed. Additionally, while meshes represent only the surface of many neighboring voxels, they do so by using three floating point numbers for the vertices plus three integers for each triangle face. For some geometries, this can inflate the size of the representation compared with the image, especially when the image is well-compressed. Skeleton fragments and final files (both 300 MB) are lightweight as expected.

6.8 GB of final compressed files remained after processing. An additional 4.6 GB of intermediate mesh fragment files and 300 MB of intermediate skeleton fragment files were also present which can be deleted. If no intermediate files are deleted, the total size of the finished dataset is 12 GB including the spatial index (which can be useful for end-users beyond the generation of meshes and skeletons).

To help illustrate Igneous's ability to scale, we provide some illustrative examples of the performance of several large jobs that have been run over the years. These jobs were run on Google Cloud Platform and Google Cloud Storage using preemptible 16-core n1-standard-16 (104 GB RAM), n1-highmem-16 (128 GB RAM), and e2-highmem-16 (128 GB RAM) machine types with one Igneous process per a core. Igneous's performance may have improved since these jobs were run.

Note that while we provide core-hours here, the total cost of these jobs depended on many factors such as the number of reads, writes, and bandwidth. Often the most expensive cost was writing unsharded meshes or skeletons as this would generate at least one or two files per a segmentation label resulting in hundreds of millions or billions of files generated. Hence, for large datasets we now recommend the sharded format.

Table 3 was compiled retrospectively from contemporaneous notes. It shows that Igneous scales to jobs requiring at least 16,000 cores and 1,000 machines simultaneously for grayscale image downsampling on a 95 teravoxel volume which was completed in a little over an hour in wall-clock time. The largest volume processed was 298 teravoxels which was a segmentation downsampling job and was completed in a little over 3 h wall-clock time.

The most computationally expensive job listed was the 5/4/19 unsharded primary mesh fragment generation step for an unpublished draft automatic segmentation of the cubic millimeter dataset (MICrONS Consortium et al., 2021) which took 99,800 core-hours. By comparison, the following merging step was much less intensive and took only 2,450 core-hours.

Skeletonization was similarly demanding. The most computationally expensive job listed was the 1/21/20 sharded skeleton fragment generation at an estimated 91,000 core-hours for the same dataset. The skeletonization parameters were set to capture dendritic spines. That run produced 397,770,063 skeletons in 524,288 shard files which occupied 747.6 GB of disk space. An example skeleton can be seen in Figure 9. The corresponding merging operation is not shown in Table 3 as it was run almost a year later and at the time required an iterative development cycle to increase performance over a period of months making a point estimate unhelpful.

The skeletonization run on 3/24/22 was conducted on an unpublished dataset with spine capture and “avocado fixing” enabled. This caused the primary skeletonization phase to be slower than would otherwise be expected. The merging run had the somewhat expensive “remove ticks” feature disabled, and thus preserved small extensions.

Though they are primarily for visualization, meshes may be used to derive scientific insights. Therefore, we endeavored to provide some basic characterization of mesh quality. Using the same CREMI A padded segmentation used in our image compression experiments, we compared Igneous/zmesh generated and stitched meshes to meshes generated by scikit-image⁵² using skimage.measure.marching_cubes version 0.18.1 using the lewiner algorithm. An example object that was meshed using both methods can be seen in Figure 10.

There are 37,366 unique labels and 37,828 26-connected components in the segmentation. Of these labels, 36,902 (98.76%) are smaller than 1,000 voxels while 464 (1.24%) are greater than or equal to 1,000 voxels. We generated all meshes using zmesh as both unsimplified and simplified versions from the base resolution segmentation. Simplification was performed using zmesh's built in algorithm using a triangle reduction target factor of 100 and a maximum error tolerance of 40 nanometers. scikit-image had a limitation where only objects 2 × 2 × 2 or larger could be meshed, so only 586 meshes could be generated.

We then used Trimesh⁵³ version 3.10.8 to check the resultant meshes for topological integrity using the mesh.is_volume property. According to their documentation, this checks watertightness (every edge is included in two faces), having consistent winding (each shared edge going in an opposite direction), and outward facing normals. All Igneous/zmesh meshes passed this check. scikit-image meshes failed this check. It is possible that they contained degenerate triangles as they failed a convexity test. Visually, spot checks of the unsimplified zmesh and scikit-image meshes were almost indistinguishable and overlapped almost exactly in space.

In order to quantify the degree to which the zmesh meshes fairly represented the underlying segmentation, we computed the volume and centroids of all labels and plotted histograms of the ratio of mesh volume (computed with Trimesh) to voxel volume for both simplified and unsimplified meshes as seen in Supplementary Figure 1.

In this figure, it can be seen that small label volumes are often grossly underestimated by the mesh while large labels are usually underestimated within 5 or 10% of the voxel volume for unsimplified and simplified meshes respectively. Marching Cubes cuts voxel corners to create a reasonable manifold, so it is sensible that small meshes will show larger deviations while larger meshes will show smaller deviations. In the bottom right of the figure, it can be seen that many very small objects get simplified to near or actually zero volume.

In Supplementary Figure 2, we computed the difference in voxels between mesh centroids and image label centroids for zmesh unsimplified and simplified and scikit-image meshes for all labels that were valid for scikit-image. We then evaluated the labels using connected-components-3d (Silversmith, 2021) and the Trimesh centroid method (which does not rely on correct manifold topology). The maximum error for zmesh was 53.1 voxels, while for scikit-image it was 417.9 voxels. The mean error for zmesh is 4.9 with a standard deviation of 6.7 for both unsimplified and simplified meshes. For scikit-image, the mean is 8.4 and the standard deviation is 27.3. In the figure, it can be seen that all three groups are fairly similar, but scikit-image has a long tail of large errors.

Characterizing skeleton quality is somewhat more difficult than with meshes due to the different skeletonizations that can be proposed for a given object depending on the needs of the user. Therefore, skeletons are somewhat more subjective though there are proposed definitions for canonical skeletons based on a grass-fire analogy, the centers of maximally inscribed spheres, and other representations (Tagliasacchi et al., 2016).

We attempted to characterize Igneous produced skeletons by coarsely comparing them with other automatically traced skeletons on a well-known dataset. Unfortunately, datasets with manually traced skeletons generally do not have a per-voxel segmentation and vice-versa. Therefore, we decided to compare automatic skeletonizations of CREMI A as was done for meshes. We attempted to locate a TEASAR implementation that was able to be installed on our available hardware, that we were able to operate properly, and was not written by one of this article's authors. However, we were unable to do so. Therefore, we made our comparisons to skeletons generated by the popular binary image skeletonization procedure in Fiji (Schindelin et al., 2012), which implements the 1994 voxel thinning algorithm by Lee and Kashyap (1994).

We used Igneous to process CREMI A at 16 × 16 × 40 nm³ resolution using the parameters const 50, scale 3, soma-accept 3500, soma-detect 1100, soma-const 300, and soma-scale 1 on all segments with greater than or equal to 1,000 voxels. We did not use the short extension (“tick”) elimination feature (though it may have slightly improved some skeletons). For Fiji, we processed the segmentation using the same size threshold at the same resolution into a series of binary TIFF files and then batch processed them. The resultant thinned binary images were processed into SWC files and then converted into Neuroglancer Precomputed skeletons that were correctly offset into the same space as the Igneous skeletons. We visually confirmed via spot checks that both sets of skeletons appeared in Neuroglancer and seemed on-balance reasonable in their topology and location in space.

We then made several comparisons between these skeletons to characterize them. First, we used the Trimesh 3.10.8 library to check the number of points that lay outside of their enclosing (zmesh unsimplified) mesh. Nine skeletons were unable to be compared due to Trimesh repeatedly crashing during the computation. In total, 332 segments were able to be compared out of 341. 0.3% of all vertices for both Fiji and Igneous skeletons lay outside the mesh. The existence of this small quantity may be due to small differences between the 16 × 16 × 40 nm³ resolution and the meshes created at 4 × 4 × 40 nm³.

In Supplementary Figure 3, we compared the difference in centroids, ratio of cable lengths, and difference in number of terminal points (vertices with only one connecting edge) for each set of skeletons. It can be seen that the maximum difference between centroids is 2,172 nm, though most are much less than that. The average distance between centroids is 188 nm with a standard deviation of 244 nm. We visually inspected the twelve segments more than 676 nm (two standard deviations larger) than the mean to determine their issues. In four cases, the thinning algorithm only skeletonized one of multiple connected components, leading to a much shorter cable length. In six cases, the thinning algorithm created a complex structure we referred to as a “beehive” (see Supplementary Figure 4) that added extraneous path length. Three cases were more ambiguous as to which was the better skeleton, but the thinning algorithm preserved more branches and holes. A more typical case where both thinning and Igneous created reasonable skeletons can be seen in Figure 11.

The mean cable length was 4778.7 nm (stddev 3358.5 nm) for Igneous skeletons and 6613.5 nm (stddev 5584.4 nm) for thinned skeletons. 281 (85.49%) Igneous skeletons were shorter than their thinned skeleton counterpart and 51 (15.41%) were larger. The outlier on the right hand side of the middle subplot of Supplementary Figure 3 is a near spherical object that was skeletonized across the diameter by Igneous but was reduced to almost a point by thinning.

As can be seen in the last panel, the voxel thinning skeletons had many more terminal points. A terminal point is a vertex that has only one edge connecting it to the rest of the skeleton. In total, the thinned skeletons had 3541 terminal points while the Igneous skeletons had 1169 (or a 303% difference). The average Igneous skeleton had 3.46 terminal points (stddev 9.65), while thinned skeletons had on average 10.67 points (stddev 3.71).