TeraStitcher - a tool for fast automatic 3D-stitching of teravoxel-sized microscopy images

Abstract

Background: Further advances in modern microscopy are leading to teravoxel-sized tiled 3D images at high resolution, thus increasing the dimension of the stitching problem of at least two orders of magnitude. The existing software solutions do not seem adequate to address the additional requirements arising from these datasets, such as the minimization of memory usage and the need to process just a small portion of data.

Results: We propose a free and fully automated 3D Stitching tool designed to match the special requirements coming out of teravoxel-sized tiled microscopy images that is able to stitch them in a reasonable time even on workstations with limited resources. The tool was tested on teravoxel-sized whole mouse brain images with micrometer resolution and it was also compared with the state-of-the-art stitching tools on megavoxel-sized publicy available datasets. This comparison confirmed that the solutions we adopted are suited for stitching very large images and also perform well on datasets with different characteristics. Indeed, some of the algorithms embedded in other stitching tools could be easily integrated in our framework if they turned out to be more effective on other classes of images. To this purpose, we designed a software architecture which separates the strategies that use efficiently memory resources from the algorithms which may depend on the characteristics of the acquired images.

Conclusions: TeraStitcher is a free tool that enables the stitching of Teravoxel-sized tiled microscopy images even on workstations with relatively limited resources of memory (<8 GB) and processing power. It exploits the knowledge of approximate tile positions and uses ad-hoc strategies and algorithms designed for such very large datasets. The produced images can be saved into a multiresolution representation to be efficiently retrieved and processed. We provide TeraStitcher both as standalone application and as plugin of the free software Vaa3D.

Figure 1

The stitching pipeline. The central column lists both automatic and optional manual steps that compose the pipeline. For each automatic step are shown alongside the corresponding command line arguments as well as their values used in a sample workflow. The left column shows the metadata flow and the most significant portions of XML files for each step. On the right column there are some figures illustrating the elaboration of both data and metadata and the number of data readings and writings. The figure which illustrates the manual refinement after the Displacements thresholding step was obtained by using the MATLAB scripts we provide for a more comprehensive analysis of metadata. Stitchable stacks are marked in blue (dark grey) and non stitchable ones are marked in yellow (light grey).

Figure 2

The software architecture of TeraStitcher. The architecture of the proposed TeraStitcher tool is composed by four modules organized in three layers with growing abstraction levels. At the lowest abstraction level, the IOManager and XML modules contain input/output routines for data and metadata, respectively. The middle layer contains the VolumeManager module, which is responsible to model data organization and to access both data and metadata through functionalities provided by the lower layer. At the top layer, the Stitcher module implements the stitching pipeline.

Figure 3

Class diagrams of VolumeManager and Stitcher modules. A StackedVolume contains as many Stack objects as the number of stacks that compose the volume and it is responsible to model data organization. The Stitcher methods implement the stitching pipeline and the strategies to use efficiently system resources, by accessing data through the StackedVolume. These strategies use the algorithms TilesPlacementAlgo and PairwiseDisplAlgo that may depend on the specific characteristics of images. Note that TilesPlacementAlgo uses Displacements produced by PairwiseDisplAlgo. Classes MST, MIP-NCC, and MIP-NCC-Displacement implement the actual algorithms used in our implementation.

Figure 4

Multi-MIP-NCC method. Each tile is split along the D direction into substacks of n_slicesslices. For each couple of adjacent substacks, three MIPs along the three directions are extracted from the overlapping regions. Then a 2D-NCC map is computed for each pair of homologous MIPs, so obtaining two displacements for each direction. Two reliability measures extracted from the NCC map are used to select for each direction the best of the two available displacements.

Figure 5

Optimal tiles placement using a minimum spanning tree (MST). Weighted undirected graph of tiles whose edges represent displacements. Weights are computed as the inverse of displacements reliability measures, yielding + ∞for the unreliable ones. These prevents the tree connecting stitchable tiles (marked as blue squares) from passing through the nontistchable ones (marked as yellow circles). If different reliability measures are available along the three directions, as in the case of using MIP-NCC to compute misalignments, three minimum spanning trees are obtained for each direction separately. The marked tile in the figure is a nonstitchable tile which is excluded from any MST path traversing stitchable tiles.

Figure 6

Result of stitching at the border of four overlapping tiles. Merging without aligning nor blending (a) and with aligning and blending (b).

Figure 7

Memory management in the pairwise stacks displacements computation step. The whole volume is processed a layer at the time, each composed by n_slicesslices. For each layer, stacks are read row-wise when there are more rows than columns (i-iv). A new substack is loaded (i) when its North substack needs it for displacement computation (ii). After that, the North substack can be released (iii) and this process is repeated for the next column (iv).

Figure 8

Memory management in the merging step. The whole volume is processed a layer at the time, each composed by the number of slices needed to generate the lowest resolution (e.g. 32 slices for a reduction of 5 times). For each layer, stacks are read row-wise and merging is performed slice by slice. When one slice group of every tile in row i of the tile matrix have been read, adjacent tiles in the groups are merged using the blending algorithm described in the Algorithms section. This produces a horizontal stripe, which is merged with the ones previously generated from the slices of the preceding rows, so producing larger stripes until the whole slice of the final image representation is obtained.