Possum—A Framework for Three-Dimensional Reconstruction of Brain Images from Serial Sections

Techniques based on imaging serial sections of brain tissue provide insight into brain structure and function. However, to compare or combine them with results from three dimensional imaging methods, reconstruction into a volumetric form is required. Currently, there are no tools for performing such a task in a streamlined way. Here we propose the Possum volumetric reconstruction framework which provides a selection of 2D to 3D image reconstruction routines allowing one to build workflows tailored to one’s specific requirements. The main components include routines for reconstruction with or without using external reference and solutions for typical issues encountered during the reconstruction process, such as propagation of the registration errors due to distorted sections. We validate the implementation using synthetic datasets and actual experimental imaging data derived from publicly available resources. We also evaluate efficiency of a subset of the algorithms implemented. The Possum framework is distributed under MIT license and it provides researchers with a possibility of building reconstruction workflows from existing components, without the need for low-level implementation. As a consequence, it also facilitates sharing and data exchange between researchers and laboratories.


Introduction
In modern neuroscience we use multiple imaging techniques to study brain structure and function. The usage of complementary approaches is considered indispensable in characterizing the spatial organization of neuronal structures and their circuitry (e.g. Annese 2012; Leergaard et al. 2012;Osten and Margrie 2013). Despite exciting developments in brain imaging there is still no perfect technique which would allow for a comprehensive insight into all aspects of the brain over a wide range of spatial and temporal scales. From a practical point of view this means that data from non-destructive and destructive imaging methods have to be integrated.
The non-destructive techniques are methods to evaluate the properties of a brain without causing damage in a way allowing further experimental procedures on a specimen. The most common examples are different varieties of magnetic resonance imaging (MRI) or diffusion tensor imaging (DTI). These techniques allow one to obtain inherently three-dimensional (3D), virtually undistorted images, which deliver quantitative information on macroscopic tissue properties (e.g. Johnson et al. 2012), and can be applied in vivo.
Despite their versatility, such techniques suffer from several drawbacks. The most important are lower imaging resolution and lower specificity in comparison with methods based on microscopic imaging of sections stained with various techniques or injected with tracers.
The destructive techniques of brain imaging are methods which once applied prevent most further experimental procedures. They rely on obtaining series of two-dimensional images of brain tissue while sectioning or once the sections are cut and mounted on slides. Due to their specificity and microscopic scale of imaging they facilitate studies of fine properties of brain tissue at the cellular and subcellular levels. Typical examples are description of cyto-and chemoarchitecture of the neural tissue (e.g. Zilles 1985;Hof and Sherwood 2005) or axonal-level connectivity (e.g. Rosa et al. 2009;Zaborszky et al. 2015) via injections of different tracers. Other examples include autoradiography (e.g. Hess et al. 1998;Lebenberg et al. 2011), imaging unstained sections (Palm et al. 2010;Annese et al. 2014), or imaging sections which underwent in-situ hybridization process (Jagalur et al. 2007; Morris et al. 2010).
The capabilities of the 2D sectioning techniques have been elevated in recent years by high-throughput processing of sections and imaging methods (e.g. Chung et al. 2011;Ragan et al. 2012;Osten and Margrie 2013) along with relevant computational routines. This made it possible to conduct large-scale projects on gene expression (Lein et al. 2007), connectivity (Oh et al. 2014), and relating brain changes to behavior (Vousden et al. 2014;Kim et al. 2015).
The main drawback of sectioning approach is the necessity of physical dissecting of the brain which means it looses its three-dimensional integrity (e.g. cannot be sliced in a different plane). This is an undesirable side effect, since the brain is three dimensional and the imaging data should also be analyzed in a three-dimensional context. Additionally, comparison of section-based imaging data with results of 3D imaging methods requires the former to be brought into volumetric form (e.g. Dauguet et al. 2009;Johnson et al. 2010) which calls for adequate computational methods. Therefore, the integration of stacks of 2D images of stained sections into volumetric form is an indispensable aspect of modern neuroimaging and relevant methods and software for section alignment are actively developed.
The process of reconstruction of a series of 2D images into volumetric form is considered a difficult and time consuming task due to tissue distortions introduced during processing of the experimental material (Breen et al. 2005;Dauguet et al. 2007). These typically include global shrinkage and dehydration due to fixation in formaldehyde, freezing or paraffin embedding. Moreover, during cryosectioning, mounting and staining, additional distortions are incurred. Shearing, tearing and displacement of individual parts of the sections, non-uniform shrinkage due to chemicals used during staining procedures are common artifacts. Workflows aiming to faithfully reconstruct the 3D image of the brain have to account for such distortions (Qiu et al. 2011).
In its simplest form, the reconstruction can be performed by sequentially aligning consecutive sections to a reference section, usually the middle one (Kiessling 2011, p. 329). However, this approach is known to be volatile and sensitive to alignment imperfections and sections' distortions (e.g. Nikou 2003). Moreover, the reconstructions performed this way tend to deviate from the true, anatomical shape. This phenomenon was nicknamed the banana effect (Malandain et al. 2004) and can be overcome by introducing a shape prior-a reference 3D image, which can be an MR image, collection of images of the face of the cutting block, atlas delineations, or a set of fiducial markers, which enforce the anatomical shape of the reconstruction.
Among the many reconstructions attempted, the recent work on the whole human brain stands out (Amunts et al. 2013;Annese et al. 2014). In both cases photographs of the face of the cutting block were taken while cryosectioning the brain. Individual sections were stained for cell bodies and then affinely aligned to corresponding blockface photographs which allowed to recover the anatomical brain shape. Additionally, the sections were nonlinearly corrected for section specific distortions which produced reconstruction of higher quality and made further coregistration to an MR image easier. Adler et al. (2014) conducted a reconstruction of stained sections of human hippocampus using post mortem MR images as a shape prior. The process was first performed with multiple series of affine transformations followed by a correction of spatial artifacts performed by warping sections to their immediate neighbors and to the MR image simultaneously. Similarly, Chakravarty et al. (2006) performed a reconstruction of human basal ganglia and thalamus by applying affine alignment step followed by nonlinear corrections based on simultaneous warping towards both immediate neighbors of a given section.
Similar, multi stage reconstruction strategies were applied in non-human primate projects. For instance, in Choe et al. (2011) myelin-stained sections of owl monkey brain were registered to two-dimensional blockface images using a combination of linear and nonlinear methods prior to linear and nonlinear registration of the blockface volume to the MR image. An analogous strategy was used by Dauguet et al. (2009) in order to create a 3D digital atlas of the thalamus based on a series of stained histological sections of a baboon brain.
In reconstructions of rodents' and other small animals' brains similar strategies are employed. Some reconstructions are performed with the help of blockface images to which corresponding histological sections are aligned affinely (e.g. MacKenzie-Graham et al. 2004;Bertrand and Nissanov 2008;Mailly et al. 2010). Other approaches, like alignment based on landmarks, are also used (e.g. Hess et al. 1998).
Looking through the literature it seems that the present methodology crystallized around multistage workflows involving some form of a shape prior as an intermediate modality to which, after preprocessing, images are aligned affinely. Afterwards, nonlinear corrections are applied to account for distortions of individual sections. The resulting volume is then registered to a reference template which is either histology-or MR-based. Such strategies were utilized in, for instance, Lein et al. (2007), Uberti et al. (2009), andLebenberg et al. (2011).
While it may seem that the reconstruction of images of serial sections should be a routine and standardized procedure, this is not the case. The above mentioned projects had specific goals, relied on different data modalities, used various numbers of specimens, etc. Additionally, the majority of the described reconstructions were done using workflows tailored to the goals of a particular project and appropriate software was not released publicly. Therefore, there is a lack of generic tools for handling this type of three dimensional reconstruction from serial sections. The available image registration tools for both 3D and 2D images (e.g. Klein et al. 2010;Avants et al. 2011;Peng et al. 2011) as well as for series of images (Thévenaz et al. 1998;Ribes et al. 2010;Cardona et al. 2012;Wang et al. 2015) perform excellently in their respective domains but they are not sufficient for reconstructions based of the whole brain histology for higher animals, from rodents to human and non-human primates, in a way mentioned in the earlier part of the introduction. A good example of these difficulties can be seen e.g. in Fig. 5 of Wang et al. (2014). The software proposed therein as a solution, despite implementing several useful methods (Wang et al. 2015), does not incorporate a reference image in the reconstruction process. Therefore, the issue of reconstructing section cut in an arbitrary plane (e.g. oblique) cannot be addressed and the banana effect cannot be mitigated.
In this article we present the Possum volumetric reconstruction framework, which is open software providing building blocks for constructing computational pipelines for reconstruction of 3D images based on serial sections. The software was created by selecting frequently utilized components of such workflows including: 1) naive affine sequential reconstructions, 2) affine sequential reconstructions designed to counteract the propagation of the artifacts, 3) routines for reconstruction with the presence of a shape prior, 4) deformable reconstruction workflow intended to account for distortions of individual sections. All workflows were implemented assuming that slices are ordered and have a constant thickness.
In addition to implementing the routines, we also validate the software using synthetic datasets illustrating properties of the provided algorithms. Additionally, we test the framework against publicly available MR and histologybased datasets. The validation results are reproducible and available as a part of the framework.
The problem addressed by the proposed software is the reconstruction of 3D images from series of stained sections. There is no attempt to handle coregistration of multiple 3D images as such a task is successfully addressed by existing software for intensity-based (e.g. Klein et al. 2010;Avants et al. 2011) or landmark-based (e.g. Peng et al. 2011) 3D image registration software.

Methods
The framework has been developed using multiple technologies and comprises a collection of fundamental workflows used to generate 3D reconstructions and accompanying image processing routines. The body of the framework was implemented in the Python programming language (http://www.python.org) while the examples are available as shell scripts for easy customization and interaction with Linux-based operating systems.
The Possum framework relies on the Advanced Normalization Tools (ANTS, RRID:nlx 75959, Avants et al. 2011) software for conducting image registration. The ANTS package is used to perform both affine (including rigid) and deformable types of registration. The latter are carried out with the symmetric image normalization (SyN) method of Avants et al. (2008). Three image similarity metrics are used in the framework: crosscorrelation coefficient (CC, Avants et al. 2008), Mattes' mutual information (MI, Mattes et al. 2001) and mean square intensity difference (MSQ, Schroeder 2005) defined as: where A i , B i are the intensities of the i-th pixels of images A and B, and N is the number of pixels comprising the images. The Neuroimaging Informatics Technology Initiative file format (NIfTI, http://nifti.nimh.nih.gov/nifti-1) was selected as the data format of choice. It is capable of holding spatial information about the image (origin, spacing, orientation, anatomical directions, etc.) and storing 2D and 3D images using multiple data types (e.g. 8 bit, 16 bit, double and 24 bit per pixel RGB images which are used in the framework).
All command line software is invoked from Python scripts using a set of wrappers which allows for trivial parallelization using GNU Parallel (Tange 2011). Workflows were designed for execution in the parallel mode on either multi-core machines or on a computational clusters under supervision of a resource manager like TORQUE (http://www.adaptivecomputing. com/products/open-source/torque/) or SLURM (https:// computing.llnl.gov/linux/slurm/).

Graph-Based Affine Sequential Alignment
Graph-based sequential alignment (Yushkevich et al. 2006a) is a procedure for 3D reconstruction intended to minimize accumulation of registration errors and to identify highly distorted sections which might disturb the reconstruction process.
The implementation provided in the Possum framework is presented in Algorithm 1. In the first step every image is affinely aligned to neighboring sections towards either end of the stack. The CC metric with values rescaled to −1, 0 (the lower the value the more similar are the images) between each coregistered pair is recorded. Afterwards, a weighted graph (G) is built with vertices (V) representing individual images and edges (E) representing weights based on similarity between the images.
A reference section I r is then selected, usually from the middle of the stack, and a transformation from I r to any other image I i is obtained by computing the least cost path W r,i in the graph G using the Dijkstra algorithm (https:// networkx.github.io/). This corresponds to a chain T r,i of partial affine transformations t i,j to be composed.
The number of chained transformations might be shorter than the nominal distance between the sections I r and I i . This is interpreted as skipping those sections which are difficult to align to their neighbors. The preference to skip outstanding sections is adjusted with parameter λ which modulates the edges' weights. Small positive λ favors section skipping while larger tends to preserve sections from being omitted in the transformation chain. Note that naive sequential alignment is a special case of this workflow for = 1 regardless of the λ value (Kiessling 2011).

Iterative Affine Pairwise Alignment
The iterative affine pairwise alignment workflow allows one to construct a 3D image from a series of serial sections in the presence of a shape prior. It is particularly suitable for reconstructions in which the cutting plane of sections does not match the corresponding plane in the reference image, e.g. when coronal sections were cut at an angle with respect to the coronal plane defined in an atlas (Malandain et al. 2004;Yang et al. 2012;Adler et al. 2014).
This procedure simultaneously improves alignment of the reconstruction to the reference image by calculating a global 3D affine transformation between the reconstruction and the reference. At the same time it aligns the experimental sections to the corresponding virtual reference cuts (see Algorithm 2). Eventually, the algorithm converges to a reconstruction which is affinely aligned to the the reference image and in which the sections are broadly consistent between one another.

Coarse-to-fine Alignment
The coarse-to-fine reconstruction approach proposed by Yushkevich et al. (2006a) uses a reference image (shape prior) in which each experimental section has a corresponding section from the reference image assigned. The method is intended to account for the accumulation of registration errors, a Z-shift (the banana effect), and to mitigate severe discontinuities in the reconstruction.
The overall idea of the approach is to perform and then combine two series of rigid registrations to produce a faithful reconstruction. The first, coarse-scale registration, relies on aligning images being reconstructed to corresponding sections of the reference image. Such sections can be obtained, for instance, by using the Iterative Affine Pairwise Alignment and then resampling the transformed reference image in the space defined by the stack of histological images. This series of transformations recovers the overall shape of the brain but does not yield accurate section-to-section alignment.
The second series of transformations, the fine-scale registration, is realized by any kind of sequential alignment workflow (e.g. naive or graph-based), and aims to provide pairs of neighboring sections well aligned to one another. A Z-shift might be introduced in this stage and the overall shape of the reconstruction might be different from the reference one.
To obtain the final result, the high-frequency component of the fine-scale alignment is combined with the coarse-scale registration. This is done by Gaussian smoothing of individual parameters of the fine-scale transformation (translation and rotation angle) across the z (stack) dimension and filtering them out before combining with the parameters of the coarse-scale registration. This yields a reconstruction which preserves both the global shape of the brain and local anatomical details, combining advantages of both coarse and fine registration.

Deformable Reconstruction
The method for deformable refinement of the reconstruction of the histological volume stems from an assumption that a change of shape of a brain structure is slower than the section thickness. Thus the neighboring images are similar to one another in a formal sense (Chakravarty et al. 2006;Ju et al. 2006;Adler et al. 2014). Analysis of the theoretical properties of the method was provided by Gaffling et al. (2014).
The elementary step of the method consists of registration of the image of a given section M i to a fixed image F i obtained by averaging images of sections in -neighborhood of section M i . Such an alignment is performed for each image in the stack which constitutes a single iteration (Algorithm 3). An arbitrary number of iterations can be conducted and the deformable registration parameters may vary between consecutive iterations.
This process eliminates high frequency discontinuities caused by section-specific distortions effectively separating the anatomy from the deformation (Gaffling et al. 2014). This translates into smoother and easier to distinguish anatomical structures.

Coarse-to-fine Workflow
To illustrate the properties of the coarse-to-fine reconstruction workflow, a training T2-weighted MR image with an isotropic voxel resolution of 1 mm 3 of a curved banana was used (I S , Fig. 1a). The image was distorted by randomly translating and rotating each of the 200 individual banana slices (Fig. 1b). The amount of translation was drawn from a Gaussian distribution with μ = 0 mm, σ = 10 mm for translation in each direction and μ = 0 • and σ = 10 • for rotation. Then, a coarse-to-fine reconstruction was performed during which the distorted image was reconstructed using the undistorted image I S as a shape prior. In the coarse-scale step, the images were aligned rigidly to the corresponding section of the undistorted image, and MI similarity metric was used. The fine-scale transformation series was calculated using the naive sequential alignment, with 110 th section designated as the reference, rigid alignment and MI similarity metric. The transformation merge was performed with σ = 5 sections for both translation and rotation parameters of the rigid transformations.
The results of the coarse-scale reconstruction (I C , Fig. 1c) show that it recovered the true shape of the phantom, however, the neighboring sections are only roughly aligned to one another; MSQ(I S , I C ) = 55. The fine-scale step (I F , Fig. 1d) resulted in reconstruction in which sections were well aligned to one another although with notable z-shift and volume twist; MSQ(I S , I C ) = 663.
The merge of the two transformation series (I M , Fig. 1e) provided a reconstruction which preserves the true  (Fig. 1f).

Graph-Based Affine Reconstruction
An MR image of the naive-sequentially aligned banana (Fig. 1d) was used to prepare a synthetic dataset to assess the efficiency of the graph-based sequential alignment. The image was distorted by applying randomly generated affine transformations independently to each section (Fig. 2a). Additional distortions were introduced by manually removing approximately a half of the banana slice on 15 sections. In particular, groups of three and then two successive distorted sections were created in this way.
The distorted image underwent the graph-based sequential alignment with various reconstruction settings. Tested values of the neighborhood radius were = 1, 2, 5 and λ values of 0 and 0.5 were used (Fig. 2b).
The results of the reconstruction (Fig. 2c, d) depend on values of both parameters, and λ. The primary difference is the amount of the recovered shape for different values of the parameter. For = 1, which is equivalent to the naive sequential alignment, highly distorted sections cannot be omitted regardless of the lambda value (see = 1 in Fig. 2c, d). This is enhanced when several heavily distorted sections follow in a row. Increasing the neighbor size to 2 makes it possible to handle two successive distorted sections. Consequently, further increment of the value allows the method to detect and omit more distorted sections for both tested λ values ( = 2, 3 in Fig. 2c, d). Comparing the reconstructions performed with different λ values confirms that the higher the λ the fewer distorted sections are skipped which is in accordance with the assumptions of the method and curves shown in Fig. 2b.

Deformable reconstruction
The dataset used to illustrate the deformable reconstruction routine was a T2*-weighted MR image of an 80 days old Wistar rat brain (Johnson et al. 2012). The original image was downsampled to 25 × 50 × 50 μm 3 voxel size which corresponds to 1600 × 400 × 400 (coronal, sagittal, horizontal planes, accordingly) voxels. The image was sliced in the coronal plane and the in-plane resolution of 50 × 50 μm 2 was preserved while the thickness of the synthetic coronal sections, d, ranged from 20 to 100 μm in the intervals of 5 μm in different reconstruction trials. The reference images obtained with this procedure for a given section thickness d will further be denoted by R d .
Subsequently, the synthetic coronal sections were nonlinearly distorted to mimic deformations naturally occurring during the preparation of the histological sections. The distortions were modeled by an application of a 2D displacement field T σ , with σ characterizing spatial correlations. Each component of the displacement field was constructed from a 2D white noise image smoothed with a Gaussian filter with kernel size of σ = 300 μm rescaled to −r, r . The deformation amplitude r was chosen so as to set the median magnitude of the displacement vector T σ = 50 μm which is a value obtained for Nissl-stained sections (Majka 2014, p. 47). The procedure was repeated for all sections in the stack yielding a 3D image with individually distorted coronal sections denoted by D d . Such an image underwent the procedure of deformable reconstruction.
We used the ANTS software to perform registration between individual images in the deformable reconstruction process using the following procedures: SyN transformation model with the gradient step of 0.025 and CC similarity metric with a kernel size of 2 voxels; Gaussian regularization with a sigma of 1 voxel for both similarity gradient and displacement field; five level multi-resolution registration scheme with 1000 iterations at each level. These settings remained unchanged for all trials. Each reconstruction trial consisted of 20 iterations. The neighborhood parameters were = 1, 2, . . . , 10, and the synthetic coronal section thickness was d = 20, 25 . . . , 100 μm. For each pair of parameters ( , d) three reconstruction trials were conducted which amounted to 510 trials. Reconstructed image after iteration i corresponding to distorted image D d is denoted as B d, ,i .
Next, we studied the reconstruction accuracy. To determine how close the reconstructions were to the initial image, for each reconstruction trial we computed the MSQ similarity measure between the reference image R d , the distorted image D d , and the reconstructed image B d, ,i .
In the first step we determined the values of for which the reconstructions were most similar to the original. For a tuple of three parameters (d, , i) the value of was chosen so that * = arg min MSQ (R d , B d, ,i

)| d,i .
In 46 cases out of 51 the most accurate reconstructions were achieved for * = 1, in the remaining 5 cases it was * = 2. We set = 1 for further analyses.
In the next step the number of iterations yielding optimal reconstruction was identified by choosing the iteration index i * leading to the highest improvement in comparison with the initial distorted image: (1) We will refer to the measure used here as relative similarity.
The results show that the thinner the section is (the larger the number of sections) the larger number of iterations is required to reach the most accurate reconstruction for the given section thickness (Fig. 3a). For the 20 μm sections the optimal number of iterations was 19 and it systematically lowered with increasing thickness so that for the section thickness of 100 μm (400 sections) only 5 iterations were required to get the reconstruction most similar to the undistorted image. Note that when the number of iterations grows, the reconstruction diverges from the reference image, thus one needs additional measures to identify the optimal reconstruction. In practical applications, when the undistorted image is unknown, the stop-point can be determined, for instance, either by letting a neuroanatomist to decide which reconstruction is the most satisfactory or by tracking the changes (e.g. calculating the MSQ value) between consecutive reconstructions and terminating the process once the difference reaches an arbitrarily defined value or the first local minimum.
Additionally, the reconstruction accuracy lowered as the section thickness increased. This is expressed by the values of the relative similarity which grew linearly with the section thickness (Fig. 3b). For the initial 20 μm thickness it was 0.27 and increased up to 0.55 for the 100 μm sections with the factor of 4 · 10 −4 per μm.
To assess how the deformable reconstruction workflow scales up with the different number of sections being reconstructed, the total CPU time elapsed on each reconstruction trial (consisting of twenty iterations) was recorded (Fig. 3c). The time increased from 21.2 CPU hours (1.1 hour per iteration) for 400 sections (section thickness of 100 μm) to 118 h (5.9 h per iteration) for 2000 sections (20 μm thick). On average, the time required to conduct a single reconstruction trial increased with a factor of 0.0031 CPU hours (11.36 s) per section per iteration. The tests were performed under Ubuntu 10.04 operating system deployed on a dual Intel® Xeon® E5620 (16 × 2.40 GHz logical processors) server equipped with 32 GB of RAM.
An example reconstruction conducted for the section thickness of 50 μm, = 1 and i * = 9 (yellow point in Fig. 3a) was selected for presentation in Fig. 4. We can see that the reconstruction was the most accurate in regions where the brain structure changed the least and for structures of low curvature, e.g. neocortex, thalamus, hypothalamus, midbrain, pons (empty triangles in Fig. 4e). On the other hand, structures with complex internal details and substantial curvature, e.g. cerebellar cortex, striatum, suffered from some reconstruction artifacts like segmentwise overstraightening and general difficulty in recovering fine details of curvature (filled triangles in Fig. 4e). Some of the artifacts are noticeable also in the olfactory bulb.  showing high discrepancy close to the brain outline, inside the cerebellum, olfactory bulb, and striatum, d a reconstructed image B d=50 μm, =1,i * =9 , e the difference between images presented on panels a and d showing that the discrepancy was significantly reduced after conducting the reconstruction process. Empty triangles denote regions in which the reconstruction was successful while filled triangles indicate worse performance

Evaluation on an Open Histological Dataset
To demonstrate the capabilities of the framework we performed a reconstruction based on the Waxholm Space Mouse Brain Atlas (Johnson et al. 2010). The dataset included (21.5 μm) 3 isotropic T2*-weighted MR image of the 80 days old CJ57BL/6 mouse brain and a series of 312 images of horizontal Nissl-stained sections of the same brain. To streamline the calculations and make the data convenient to share as an example, the original MR image c b a Fig. 5 Results of a 3D reconstruction of the Nissl-stained sections (Johnson et al. 2010). Column a: the outline of the 3D reconstruction of the brain (gray model) and the location of the cuboid shown in the column (b). Coronal (a1), lateral (a2), and horizontal (a3) projections. The plane shown in panel (a1) was used to obtain oblique cuts through the reconstruction and the reference image presented in column (c). Column b: Fragments of the reconstructions showing dentate gyrus and neighboring structures during consecutive stages of processing, 1) pairwise alignment, 2) graph-based sequential alignment, 3) deformable reconstruction, the final stage, 4) the reference MR image. Column c: Oblique cuts through the reconstructions: 1) the final reconstruction, 2) the reference MR image. Red arrows on panel c1 indicate successful, detailed reconstructions of olfactory bulb, piriform cortex, hippocampus, and cerebellum was downsampled to isotropic (43 μm) 3 resolution and the high-resolution Nissl-stained sections were downsampled to pixel size of 50 × 50 μm 2 . A reconstruction workflow comprising the following four steps was applied.
To begin with, ten iterations of the pairwise registration workflow were conducted. During this step the histological images were rigidly aligned to the corresponding virtual sections obtained by affinely aligning the reference image to the stack of histological images being reconstructed. Correlation coefficient (CC) was used as the image similarity metric. This stage resulted in a rough registration of the histological image stack to the reference MR image.
Subsequently, graph-based sequential alignment with = 5 and λ = 0 was applied to the results of the previous step. The image of the 110 th section was used as the reference to which all the remaining section were aligned sequentially using rigid transformations and MI as the image similarity metric.
The last step was to apply eight iterations of the deformable reconstruction workflow. During this process the following parameters were used: = 1; CC image similarity metric with the kernel size of 2 voxels; gradient step of 0.01; Gaussian regularization with kernel size of 100 μm for similarity gradient and 50 μm for the displacement field; six-level image pyramid with 1000 iterations per level.
Ultimately, transformations from all the intermediate steps for every section were merged and used to reslice this section. The reconstruction yielded 3D brain image coregistered affinely to the reference MR image (Fig. 5). Additionally, the same transformations were applied to the masks of sections which provided a reconstruction of the brain outline (Fig. 5a).

Discussion and Summary
In this article we introduced the Possum framework, software addressing the task of reconstruction of threedimensional images of the brain based on series of twodimensional images of stained sections. To develop the framework we reviewed and selected workflows which accomplish versatile reconstruction tasks according to today's best practices. Additionally, we demonstrated properties of individual routines and illustrated which reconstruction task they are suitable for.
The framework design paradigm was to utilize reliable open source components and follow guidelines for designing and evaluating scientific software (e.g. Tustison et al. 2013). This decision resulted in both increase of the framework's stability and reduction of the development efforts. The primary external component is the ANTS software, chosen because it is a thoroughly tested (Klein et al. 2009), customizable, and task-agnostic image registration tool. By using open source and well maintained image processing libraries (e.g. ITK) and NIfTI file format which are de facto standards (Poline et al. 2012;Avants et al. 2015) in neuroimaging we increased the interoperability between the framework and other pre-and post-processing tools. However, the chosen paradigm caused an overhead due to data exchange between individual components and resulted in hampering, to some extent, the efficiency of the framework. Replacing external dependencies with dedicated components is a part of undergoing maintenance work. The scalability of the software is mainly a result of enabling parallel processing without which handling sizable datasets would be inconvenient. Due to the nature of the data-multiple images which usually can be processed independentlynaive parallelization turned out to be a sufficient solution.
To assess the scalability we tested the framework against both, relatively small synthetic datasets of the banana as well as large ones, containing up to two thousand sections. The framework managed to handle both situations well.
The implemented workflows constitute a software collection which allows one to conduct typical tasks of reconstructing 3D brain images from series of stained sections assuming they are properly ordered and of constant thickness. Graph-based affine sequential alignment and coarseto-fine routines allow one to conduct affine reconstructions with or without using a shape prior while at the same time reducing Z-shift, skewing, banana effect, and propagation of the alignment errors due to the sections' distortions. The iterative alignment workflow makes it possible to perform the reconstruction within the coordinate space of the reference image accounting for the fact that the sections were cut in a different plane than the corresponding sections in the reference image. The deformable reconstruction workflow uses nonlinear transformations to compensate for sectionspecific distortions, eliminating discontinuities, improving overall reconstruction quality, and making it easier to conduct further 3D to 3D coregistration tasks. Individual reconstruction routines can be stacked to create pipelines tailored to specific projects.
With the example of Nissl stained sections of the Waxholm Space Mouse Brain we showed that the framework is capable of tackling data from demanding research projects. The workflow used there addresses a typical reconstruction task in which a 3D brain image is reconstructed from a series of stained sections and a reference image. Additionally, during its development, the framework was used to create stereotaxic atlas of the Monodelphis opossum brain (Majka et al. 2013) and to construct a workflow for connectivity data mapping in the common marmoset brain .
One aspect of 3D reconstruction not addressed by methods presented in this article is repairing highly distorted sections (e.g. detached parts of the tissue, tears, etc.). Different approaches have been utilized to mitigate such distortions. For instance, Choe et al. (2011) corrected displaced parts of sections by identifying them manually in, both, histological and reference (blockface) images. A more automated approach was proposed by Dauguet et al. (2007) where semi-automatic hemisphere separation was performed assuming that the processed sections were symmetric. A more elaborate approach was used by Amunts et al. (2013), who iteratively perform 3D reconstruction and remove minor artifacts while severe distortions are still identified and corrected manually.
For detached and torn pieces of tissue methods similar to those proposed by Pitiot et al. (2006) and Pitiot and Guimond (2008) seem to be a good remedy. Briefly, their algorithms model such distortions with several rigid or affine local transformations embedded in an elastic one. Not only it allows to automate repairing heavily distorted sections but also makes it possible to encode such corrections as displacement fields which is important from the reproducibility standpoint.

Framework in the Context of Digital Brain Atlasing
The presented software addresses the key issue of creating 3D images from serial sections before deformable mapping to a three-dimensional reference space using existing software. By reducing the efforts of establishing 3D reconstruction pipeline and providing reliable routines, the Possum framework gives the researchers an opportunity to conduct projects involving histological data integration by themselves. The framework might be used for instance to facilitate delineation of brain structures based on both, MR images and histology (e.g. Kumazawa-Manita et al. 2013;Ullmann et al. 2015), or in brain connectivity studies (Kuan et al. 2015;Sukhinin et al. 2015).
Another example might be processing legacy data by which we understand histological experimental material which has been collected without intention to reconstruct in it 3D but which may still constitute a valuable neuroscientific resource and therefore would benefit from integration with other digital atlasing resources, such as Zakiewicz et al. (2015).

Further Directions and Outlook
The directions of further development are twofold. In terms of the framework functionality, the next step is to provide routines for preprocessing of the input data, e.g. interfaces for managing collections of input images, implementing routines for correcting section staining inhomogeneities (e.g. Chakravarty et al. 2006;Yelnik et al. 2007;Ceritoglu et al. 2010), fixing tears or severe displacements of the tissue (e.g. Dauguet et al. 2007;Pitiot and Guimond 2008) before conducting the actual reconstruction process. With regard to the framework architecture, the next step is to develop a mechanism for easy connecting, interfacing between consecutive steps and monitoring execution of the pipeline (e.g. Gorgolewski et al. 2011;Friedel et al. 2014). We ultimately envision the framework as a back-end of a web service connected to high resolution image repositories (e.g. Mikula et al. 2007) which would allow one to compose and execute various reconstruction pipelines.

Information Sharing Statement
The source code of the framework is distributed under the terms of the MIT license and is available to download from the GitHub repository: https://github.com/pmajka/poSSum. The repository contains also the code and data necessary to reproduce examples shown in this article.
A preconfigured virtual machine with ready-to-use installation of the framework as well as the code used to perform the reconstruction of the Waxholm Space Mouse Brain Reference can be found at: http://www.3dbar.org/ wiki/barPosSupp.