1 Introduction
Electron microscopy (EM) images contain anatomical information relevant for detailed reconstruction of neuronal circuits. Modern EM acquisition systems produce increasingly large sample volumes, with millions of images
(Zheng et al., 2017). Due to the limited fieldofview of current imaging hardware, and the desire to image large volumes, partially overlapping images (tiles) are acquired, to be registered into a seamless threedimensional volume in a postacquisition processing step. Analysis is performed on such fully registered image volumes (Takemura et al., 2015; ya Takemura et al., 2013; Bock et al., 2011; Saalfeld et al., 2010, 2012; Scheffer et al., 2013; Cardona et al., 2012).A major challenge for EM volume registration is the large dataset size. For example, to image a whole fruitfly brain (approximately 500 µm in its largest dimension) at voxel dimensions of 4×4×40 nm with serial section transmission EM (ssTEM), with 10% overlap of adjacent tiles, approximately 21 million images (380 tera voxels) are generated (Zheng et al., 2017). There is a need for robust, efficient and scalable methods for registering such volumes.
This manuscript addresses three stages of EM volume registration:
 Montage registration

The process by which partially overlapping images within the same section are registered to form a seamless mosaic (Figure 1b),
 Rough alignment

Scaleddown snapshots of whole montages produced by the previous step are roughly registered to each other across , and finally,
 Fine alignment

In which pointmatches obtained from step 1, together with pointmatches across sections are used to jointly register the whole volume.
Each of the above steps is itself a twostep process. First: Matching pointpairs are found between pairs of images identified as potential neighbors (Figure 1a and b (left panel)), based on matching image features (Lowe, 2004; Saalfeld et al., 2010, 2012). Second: Using that set of point matches, and after deciding on a transformation model, image transformation parameters are estimated by a solver algorithm. The algorithm minimizes the sum of squared distances between all point matches (Figure 1). Deforming tiles according to these estimated transformations ideally results in a seamless registration of the whole volume (Figure 1e).
Available solver algorithms for EM use pointmatch information local to each tile to find a local solution, move a step towards that solution, iterate similarly over all tiles and repeat the process until a convergence criterion is satisfied. The solver by Karsh (2016), although scalable to large problem sizes, is based on the local iterative scheme described above, regularizes against stage coordinates which are often unreliable, and is not guaranteed to converge to a global optimum. The strategy by Saalfeld et al. (2010, 2012) is limited to smaller volumes than described in this work. Restriction to iterative techniques necessitates additional solver parameterization (maximum step size, number of iterations and other empirical termination criteria) and, at the leastsquares level, limits accuracy. Other techniques, for example that by Tasdizen et al. (2010), compose pairwise transformations sequentially and are therefore prone to additive error.
Main contributions of this work:
(a) We provide a formulation of modelbased registration in linear algebra terms, to perform joint optimization of all images with explicit regularization. We thereby delegate both scalability and convergence of the solver to efficient established linear algebra computer programs that may or may not be iterative (and are generally more efficient even when iterative). For small to moderate problem sizes, our method thus uses reliable fast direct methods, which are parameter free at the leastsquares level, and outperform all iterative approaches. For large problem sizes, known iterative techniques can be used. Moreover, nonlinear models, such as higher order polynomials, and arbitrary spatial dimensions are permitted in our matrix framework.
(b) We show how, in the absence of reliable starting guesses for tile transformations, a rigid model approximation can be estimated and used for explicit regularization. The presented method therefore assumes no prior knowledge of image orientation or relative scale, but can take full advantage of priors when provided.
In practice, tiles must be corrected for nonaffine lensdistortion prior to montaging for which we used the method by Kaynig et al. (2010) followed by affine normalization. We further converted the final solution of the method described in this manuscript to a thinplate spline (Bookstein, 1989) warpfield that we applied to the seamlessly stitched montages. This way, we retain the overall shape of the deformable registration approximated by an affine transformation for each tile while avoiding discontinuous seams between tiles in the montage.
2 Methods
2.1 The leastsquares system
For a system of partially overlapping image tiles (number of tiles = ), we assume that we are given a full set of pointmatches between all neighboring tiles. We set up and solve a regularized (Tikhonovlike) leastsquares system for finding the set of transformations that optimally and jointly minimizes the square distances between all pointmatches
(1) 
The first term is the primary objective function. is a large matrix of pointpair correspondences.
is a column vector of transformation coefficients of length
, and is a column vector of length ( and are generally sparse). is a product of the total number of transformations () and the number of coefficients () for the type of transformation considered (e.g. for twodimensional affine). is the total number of pointpairs and is the individual pointmatch confidence estimate. The second term is a regularization term in which is a matrix with dimensions and is a column vector of length with values that the solution is expected to stay close to, i.e. it is an explicit regularization. only has elements on the diagonal. is a regularization parameter that controls the relative importance of fidelity of to pointmatch data vs. adherence to values in . In practice, we generalize to be a diagonal matrix which provides the ability to constrain individual parameters differently if desired.2.2 Constructing : “No regularization” case
Let us first consider solving the first term of Equation 1 without regularization, i.e. without the second term in Equation 1. First, we construct and . Matching pointpairs between overlapping tile pair are given by and , with indicating the pair (e.g.: or ) and point index . For reference see Figure 1. and are each 2vectors (with columns and ). The length of and is , and varies with the number of pointmatches found between tiles in the pair . is twice the sum of all .
For example, we want to find twodimensional transformations to register three overlapping image tiles jointly (), where each tile overlaps with the other two tiles
where and represent pointmatch matrix blocks for two connected tiles. lists point coordinates in the first tile’s coordinate system and those in the second tile’s. The structure of these blocks depends on the transformation model chosen and reflects the basis functions considered.
For example, in case of an affine model, each point in tile coordinates is transformed to the common coordinate system such that
The parameters are the transformation parameters sought. The number of columns of is , and the pointmatch blocks are given by
in which and are two coordinates corresponding to images 1 and 2, in a coordinate system local to each image. It is instructive to keep the above nonreduced system in mind, even though it does not have a solution. To solve it, one image tile (in the above case we choose tile 1) must be fixed. This sets the frame of reference for all other transformations, such that
(2) 
To fix a tile (done only for demonstration purposes here), the column range corresponding to its transformation parameters is eliminated. Vector is modified as shown in Equation 2.
For small numbers of tiles, solving Equation 2 will yield the required transformations using standard linear solvers (see Figure 1 for an example with six tiles).
2.3 Regularization
The strategy of fixing one tile is not sufficient in the general case. With increasing numbers of tiles, and when tiles are further removed from the fixed reference tile, distortions are observed at the scale of the whole layer being stitched (Figures 2 and 2c left panel). To solve this problem, we include explicit constraints, by adding a regularization term (term 2 in Equation 1).
Stagereported tile coordinates may be accurate enough to be used as an explicit constraint in special cases. However, in the general case, these numbers are not reliable. We now describe a strategy to obtain a rough guess ( in term 2 in Equation 1) to serve as regularizer. and are translated to their respective centers of mass, producing and . A two dimensional affine transformation can be constrained to a similarity deformation by considering, in addition to the pointmatch data and , the pointmatch data subjected to an operator that transforms into (Schaefer et al., 2006). We generalize this idea to the full set of pointmatches among all tiles. The system is thereby implicitly similaritydeformationconstrained by the data; both original and artificially generated. To accomplish this, we write an equation similar to Equation 2 for the joint system
(3) 
where is a column vector of coefficients (missing translation terms). The key difference is that blocks of and are now given by
When restricted to two images, this resembles the problem constructed by Schaefer et al. (2006).
We solve Equation 3 to obtain values for all transformation coefficients. Due to the requirement of fixing one tile as reference to be able to solve the matrix system, tiles that are far away from the reference tile suffer excessive reduction in scale. This is because in the limit of infinite overlapping tiles, and the fact that our deformation model is never fully accurate, error accumulates to a degree that drives the optimization to reduce this error by reducing overall scale of the tile collection.Therefore, all tiles must be subsequently rescaled to their original area to yield the desired rotation approximation . To obtain translation parameters , we solve a translationonly leastsquares system separately. The combined parameters and are used to populate column vector in Equation 1. So we write term 2 of Equation 1 as
(4) 
where has dimensions , with diagonal elements determining relative importance of regularization for a specific parameter.
is the identity matrix in most applications.
is a column vector of length and represents the approximate solution to the rigidmodel problem.Using Equations 4 and 2, and choosing a suitable value for the regularization term , we solve the full regularized system (Equation 1). In normal equation form
and assuming vector to generally be all zeros, the identity matrix and a diagonal matrix (in the general case, each parameter can be constrained independently)
so we are solving a system
(5) 
2.4 Choice of regularization parameters
For problem sizes at the scale of several thousand tiles, solution of Equation 5 is fast ( s) and a parameter sweep for determining is practical. We calculate tile transformations and plot log vs. a measure of deformation, which is taken to be the ratio of average area of deformed tile relative to area of the undeformed tile (Figure 2a).
2.5 Generalization to nonlinear transformations
The matrix expressions above extend to nonlinear models in a straightforward manner solely by modification of matrix blocks and .
2.6 Explicit constraints
Equation 5 includes the regularization parameter , which in the general case is a diagonal matrix with diagonal elements corresponding to the degree of regularization desired for each individual parameter. For example, the user might want to decrease constraints on translation parameters of each tile (leave and relatively free), and strongly constrain all parameters (including translation) for one or more sections that should not be modified by the solution. Such strategies come in handy when performing local improvements on alignment of especially problematic sections, while preventing any perturbation of neighboring sections.
3 Results
We performed a series of registrations of ssTEM image data with increasing numbers of tiles for both single section slices (montage registration) and multiplesection volumes.
In all cases, tile metadata was first ingested into a dedicated database.^{1}^{1}1Renderer: https://github.com/saalfeldlab/render Pointmatches between potential tilepairs were then calculated as by Saalfeld et al. (2012) and subsequently ingested into the database for retrieval by the solver process to build a linear system (Equation 5). Different linear solvers were used to solve Equation 5. All experiments were conducted on a 32 CPU Broadwell computer with 256 GB RAM using Matlab version 2017a. A parallel pool with all 32 CPUs was used with parfor (parallel for loops) for constructing the linear system. An explicit parallel solution in Matlab was not used. In the case of PaStiX (Hénon et al., 2002), we used a setup of 8 CPUs on a dedicated Broadwell node. Results are summarized in Tables 2 and 1. The purpose here is to provide a general idea of CPU performance obtainable for such systems with current hardware, and to compare direct vs. iterative linear solvers.
We observe that direct methods outperform iterative ones both in CPU time requirements, pertile pointmatch residual error and linear system precision as expected. For large systems with more than 1M tiles, PaStiX outperformed other approaches significantly. PaStiX is the only massively parallel linear solver that we tested. It is likely that other parallel direct solvers are also suitable for this type of problem.
The efficiency of generating montage solutions makes a regularization parametersweep practical. Equation 5 was solved for a range of values for a 536tile section dataset (Figure 2). A tiledeformation measure was determined as the mean deviation of tile areas postregistration from starting undeformed tile areas.
4 Discussion
The presented method enables joint deformable registration of millions of images using known linear algebra techniques. Leastsquares systems resembling Equation 1 are known in the literature (Golub and Loan, 2013). The main contribution of this work is enabling linear algebra solvers for the large EM registration problem by providing an explicit matrixbased formulation for joint estimation of a rigidmodel approximation. Without such a model, it is not possible to use Equation 1 for any but the most trivial problems. Importantly, parallel direct solvers can be used.
If the problem size (or hardware restriction) necessitates the use of iterative methods over direct ones, then established iterative strategies such as GMRES or stabilized biconjugate gradients may be used. In this way, image registration efforts for EM are decoupled from solver strategies and automatically benefit from existing general scalable linear solvers and future work to improve them.
4.1 Code
The computer code accompanying this work estimates transformation parameters for translation, rigid approximation, affine and higher order polynomials up to third degree. It is written in the Matlab (The Mathworks Inc.) programming language. The main solver functions are summarized in Table 3 and corresponding code can be obtained freely (Khairy, 2018).
5 Acknowledgements
We thank Davi Bock, Zhihao Zheng, Camenzind Robinson, Eric Perlman, Rick Fetter and Nirmala Iyer for generating and providing us with the full adult fly brain (FAFB) dataset. We thank Bill Karsh and Lou Scheffer for discussions about solver strategies. We thank Eric Trautman, Tom Dolafi, Philipp Hanslovsky, John Bogovic, Eric Perlman, Cristian Goina and Tom Kazimiers for building computational tools to enable data management, retrieval, zcorrection and inspection. We thank Goran Ceric, Rob Lines and Ken Carlile for help with enabling PaSTiX computations on Janelia’s compute cluster.
References
 Bock et al. (2011) D. D. Bock, W.C. A. Lee, A. M. Kerlin, M. L. Andermann, G. Hood, A. W. Wetzel, S. Yurgenson, E. R. Soucy, H. S. Kim, and R. C. Reid. Network anatomy and in vivo physiology of visual cortical neurons. Nature, 471:177–182, 2011. doi: doi:10.1038/nature09802.
 Bookstein (1989) F. L. Bookstein. Principal warps: thinplate splines and the decomposition of deformations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(6):567–585, 1989.
 Cardona et al. (2012) A. Cardona, S. Saalfeld, J. Schindelin, I. ArgandaCarreras, S. Preibisch, M. Longair, P. Tomančák, V. Hartenstein, and R. J. Douglas. Trakem2 software for neural circuit reconstruction. PLoS ONE, 7(6):e38011, 2012. doi: 10.1371/journal.pone.0038011.
 Golub and Loan (2013) G. Golub and C. V. Loan. Matrix Computations. The Johns Hopkins University Press, 2013.
 Hénon et al. (2002) P. Hénon, P. Ramet, and J. Roman. PaStiX: A HighPerformance Parallel Direct Solver for Sparse Symmetric Definite Systems. Parallel Computing, 28(2):301–321, Jan. 2002.
 Karsh (2016) B. Karsh. Aligner for large scale serial section image data. GitHub repository, https://github.com/billkarsh/Alignment_Projects, 2016.
 Kaynig et al. (2010) V. Kaynig, B. Fischer, E. Müller, and J. M. Buhmann. Fully automatic stitching and distortion correction of transmission electron microscope images. Journal of Structural Biology, 171(2):163–173, 2010.
 Khairy (2018) K. A. Khairy. EM aligner. GitHub repository, https://github.com/khaledkhairy/EM_aligner, 2018.

Lowe (2004)
D. G. Lowe.
Distinctive image features from scaleinvariant keypoints.
International Journal of Computer Vision
, 60(2):91–110, 2004.  Saalfeld et al. (2010) S. Saalfeld, A. Cardona, V. Hartenstein, and P. Tomančák. Asrigidaspossible mosaicking and serial section registration of large sstem datasets. Bioinformatics, 26(12):i57–i63, 2010. doi: 10.1093/bioinformatics/btq219.
 Saalfeld et al. (2012) S. Saalfeld, R. Fetter, A. Cardona, and P. Tomančák. Elastic volume reconstruction from series of ultrathin microscopy sections. Nature Methods, 9(7):717–720, 2012. doi: 10.1038/nmeth.2072.
 Schaefer et al. (2006) S. Schaefer, T. McPhail, and J. Warren. Image deformation using moving least squares. ACM Transactions on Graphics, 25(3):533–540, 2006. doi: http://doi.acm.org/10.1145/1141911.1141920.
 Scheffer et al. (2013) L. K. Scheffer, B. Karsh, and S. Vitaladevun. Automated alignment of imperfect em images for neural reconstruction. arXiv:1304.6034 [qbio.QM], 2013.
 Takemura et al. (2015) S.y. Takemura, C. S. Xu, Z. Lu, P. K. Rivlin, T. Parag, D. J. Olbris, S. Plaza, T. Zhao, W. T. Katz, L. Umayam, C. Weaver, H. F. Hess, J. A. Horne, J. NunezIglesias, R. Aniceto, L.A. Chang, S. Lauchie, A. Nasca, O. Ogundeyi, C. Sigmund, S. Takemura, J. Tran, C. Langille, K. Le Lacheur, S. McLin, A. Shinomiya, D. B. Chklovskii, I. A. Meinertzhagen, and L. K. Scheffer. Synaptic circuits and their variations within different columns in the visual system of drosophila. Proceedings of the National Academy of Sciences, 112(44):13711–13716, 2015. ISSN 00278424. doi: 10.1073/pnas.1509820112.
 Tasdizen et al. (2010) T. Tasdizen, P. Koshevoy, B. C. Grimm, J. R. Anderson, B. W. Jones, C. B. Watt, R. T. Whitaker, and R. E. Marc. Automatic mosaicking and volume assembly for highthroughput serialsection transmission electron microscopy. Journal of Neuroscience Methods, 193(1):132–144, 2010. doi: 10.1016/j.jneumeth.2010.08.001.
 ya Takemura et al. (2013) S. ya Takemura, A. Bharioke, Z. Lu, A. Nern, S. Vitaladevuni, P. K. Rivlin, W. T. Katz, D. J. Olbris, S. M. Plaza, P. Winston, T. Zhao, J. A. Horne, R. D. Fetter, S. Takemura, K. Blazek, L.A. Chang, O. Ogundeyi, M. A. Saunders, V. Shapiro, C. Sigmund, G. M. Rubin, L. K. Scheffer, I. A. Meinertzhagen, and D. B. Chklovskii. A visual motion detection circuit suggested by Drosophila connectomics. Nature, 500:175–181, 2013. doi: 10.1038/nature12450.
 Zheng et al. (2017) Z. Zheng, J. S. Lauritzen, E. Perlman, C. G. Robinson, M. Nichols, D. Milkie, O. Torrens, J. Price, C. B. Fisher, N. Sharifi, S. A. CalleSchuler, L. Kmecova, I. J. Ali, B. Karsh, E. T. Trautman, J. Bogovic, P. Hanslovsky, G. S. X. E. Jefferis, M. Kazhdan, K. Khairy, S. Saalfeld, R. D. Fetter, and D. D. Bock. A complete electron microscopy volume of the brain of adult drosophila melanogaster. bioRxiv, page 140905, 2017.
Comments
There are no comments yet.