Marker-free image registration of electron tomography tilt-series
© Sorzano et al; licensee BioMed Central Ltd. 2009
Received: 04 September 2008
Accepted: 27 April 2009
Published: 27 April 2009
Tilt series are commonly used in electron tomography as a means of collecting three-dimensional information from two-dimensional projections. A common problem encountered is the projection alignment prior to 3D reconstruction. Current alignment techniques usually employ gold particles or image derived markers to correctly align the images. When these markers are not present, correlation between adjacent views is used to align them. However, sequential pairwise correlation is prone to bias and the resulting alignment is not always optimal.
In this paper we introduce an algorithm to find regions of the tilt series which can be tracked within a subseries of the tilt series. These regions act as landmarks allowing the determination of the alignment parameters. We show our results with synthetic data as well as experimental cryo electron tomography.
Our algorithm is able to correctly align a single-tilt tomographic series without the help of fiducial markers thanks to the detection of thousands of small image patches that can be tracked over a short number of images in the series.
Electron tomography is a rapidly growing technique that produces structural information of organelles and cell compartments at a resolution between 40 and 20Å [1–3]. This kind of information helps structural biologists to understand how macromolecular machines interoperate in the cell to perform their functions and is giving rise to what is called "visual proteomics" [4, 5].
This technique relies on the tomographic reconstruction of the object under study from the projection images acquired by the electron microscope. The most popular data collection geometry is a single-axis tilt series , although other geometries are also possible like multiple-axis [7–9] or conical tilt [10, 11]. In this article we concentrate on single-axis tilt series. With this technique the sample is tilted around an axis to provide different projection views of the object under study. As the object is tilted electrons have to traverse thicker sections. Thus, this technique is normally restricted to a maximum tilt of about 70°. Due to mechanical instabilities of the sample holder during rotation, it is possible that some of the images of the series exhibit a large shift with respect to the rest. Besides their shift, their relative projection directions must be established before entering the three-dimensional reconstruction process. This basically amounts to determining the tilt axis in each of the images since its exact position may not be the same on each projection.
Alignment of the tilt series is usually performed by aligning gold beads that serve as fiducial markers [12–14]. Although quite effective when these markers are available, this technique is not always applicable since markers are not always visible or trackable, or because, simply, they are not available . Alternatively, image corners can be used as fiducial markers , or the image information content itself can be used through cross-correlation [17, 18]. The image information has been proven to be the most valuable information in image registration  and, therefore, algorithms based on cross-correlation with a reference volume or similar techniques are supposed to provide the best results as is the standard case in electron microscopy of single particles . However, when aligning tilt series there is no a priori volume that can be used as reference, and for this reason, the algorithm employed must align the images with information exclusively contained in the image series. There are three main possibilities to use cross-correlation in this context of tilt series alignment. The first possibility is to serially align the first image with the second by cross-correlation, the second with the third, etc. This approach compares similar images, although it has the drawback of serious potential drifts (by error propagation) of the alignment parameters resulting in inaccurate estimations of the tilt axis and shift parameters in each of the images. The second option is to use cross-correlation to track features within the tilt series. Features are defined as fiducial points in the images (due to change of contrast, borders, corners, etc.) that are searched by cross-correlation in neighboring images. The feature tracking allows to define markers that are visible in only a subsequence of the whole tilt series. The third approach constructs a rough reconstruction of the tomogram and realigns the tilt series with respect to reprojections of the rough reconstruction .
We exploit the second option by searching for image regions that can be tracked within a subseries of the tilt series. An overall affine image registration is performed to be able to predict the position of any image region in any of the adjacent images. If all images are numbered according to their tilt angle in the tilt series, a region of image i is sought in image i + 1, first predicting its position through the affine transformation, and then by local refinement of the cross-correlation. Once, the region is found in image i + 1, the equivalent region in image i + 1 is sought in image i + 2. The algorithm thus propagates a coordinate in image i to images i + 1, i + 2, until the region sought cannot be found anymore because it becomes obscured by other objects, it changes its shape, or the noise does not allow to find it. The same is done backwards. All the coordinates corresponding to equivalent image regions form a landmark chain which most likely does not cover the whole tilt series. Landmark chains are further refined to guarantee that they correspond to equivalent image regions. Finally, all landmark chains are used to perform a robust regression of the image alignment parameters. An advantage of this approach is that it is fully automatic and it can produce thousands of landmark chains which are used to align the images. We have tested our algorithm with synthetic data as well as experimental cryo electron microscopy data.
The method proposed in this paper differs from other published methods in a number of ways. The algorithm proposed by Penczek et al. in  is very similar in its basics to the one presented in this paper. Our main differences are that we provide a mechanism to automatically estimate landmark positions, our algorithm is based on the use of thousands of short landmark chains (although this feature is also allowed in Penczek's algorithm, it was not conceived to be used in this way), we perform a robust estimation of the alignment parameters, and we provide matrix closed forms of the parameter estimates which are more convenient for implementation. The algorithm described by Mastronarde in  differs also mainly in our automatic estimation of the landmark positions yielding thousands of short landmark chains. The optimization procedure is also different, ours being more analytical. The approach of Brandt in  and Castaño-Díez in  differ from ours in the way landmark chains are automatically detected and refined, and the posterior optimization procedure followed to estimate the alignment parameters. Their algorithms also perform a robust optimization, although in a different way to ours.
Results and discussion
Before reading this section we strongly suggest to read the Methods Section to get acquainted with the terminology and different parameters of the algorithm.
In order to test the accuracy of the new alignment algorithm we created a phantom built with mathematical descriptions of cylinders and spheres (see Fig. 1). The phantom size is 512 × 512 × 128 voxels. For the first experiment 41 images of size 512 × 512 were simulated with a maximum tilt of 60° and a step of 3° between images. The tilt axis passed through the coordinate origin and its direction was defined by the angles α = 0 and β = 90 (see Methods). No in-plane transformations were applied to the projections. Landmark regions had a size of 7.5% the total size of the images (regions were 39 × 39 pixels in this case), a grid of 40 × 40 evenly spaced points in each projection image were considered as centers of the landmark regions (this means that between two adjacent regions of length 39, there is an overlap of 27 pixels).
For the second experiment 61 images were simulated with a maximum tilt of 60° and a step of 2° between images. The rest of the conditions were as in the previous example except that α = 78 (this experiment was done in order to test the algorithm ability to find a tilt axis that is not aligned with any of the standard axes). The parameters for our algorithm were identical as in the first experiment. 16562 landmark chains were identified with an average length of 23.1 images. 8748 landmark chains were used in the last regression iteration. There was an average error between the projection of the 3D landmarks and its observed projections of 0.51 pixels. The correlation between the volume reconstructed with the estimated parameters and the reconstruction performed with the ground-truth parameters was 0.99.
In the third experiment, images were collected as in experiment 2 except that α = 0. The collected images were randomly shifted in X and Y by a displacement following a Gaussian distribution with zero mean and standard deviation 20 (note that the images were of size 512 × 512). 12048 landmark chains were identified with an average length of 21.38 images. 6284 landmark chains were used in the last regression iteration. There was an average error between the projection of the 3D landmarks and its observed projections of 1.13 pixels. The correlation between the volume reconstructed with the estimated parameters and the reconstruction performed with the ground-truth parameters was 0.98.
The fourth experiment was aimed at characterizing the performance of the algorithm under slowly varying drifts. The 61 images of the second experiment were used for this last experiment. The first image was rotated -4 degrees, the second image was rotated -2 degrees, the third image was not rotated, the fourth image was rotated 2 degrees, and the fifth image was rotated 4 degrees. This sequence (-4, -2, 0, 2, 4) was repeated along the rest of the series. After rotation, the first image was shifted -30 pixels in each direction. Then, the rest of the images were shifted one pixel vertically and horizontally with respect to its predecessor so that the last image was shifted 61 pixels with respect to the first one. Our algorithm was run with chain lengths of 5, a correlation threshold of 0.99, local regions of a size of 7.5% the total size of the images. 28843 landmark chains were identified with an average length of 17.74 images. 14814 landmark chains were used in the last regression iteration. There was an average error between the projection of the 3D landmarks and its observed projections of 0.33 pixels. The correlation between the volume reconstructed with the estimated parameters and the reconstruction performed with the ground-truth parameters was 0.99.
In the last experiment, we checked whether the algorithm was able to work in a noisy environment. The experimental setup was the same as in the previous experiment but we added white noise to the simulated images (the level of noise added was such that the correlation between pairs of consecutive images in the tilt series was in the same order of the correlation observed between consecutive images in the experimental tilt series described below for the Pyrodictium abyssi). The presence of noise significantly lowers the correlation between patches. For this reason we reduced the patch size to only 4% of the image size and lowered our correlation threshold to 0.9. 7848 chains were automatically extracted, out of which 4134 participated in the last regression iteration. The average of the error between the projection of the 3D landmarks and the observed projections was 2.7 pixels. Finally, the correlation between the volume reconstructed with the estimated alignment parameters and the reconstruction performed with the ground-truth parameters was 0.94.
In this article we have presented a new algorithm for the automatic alignment of electron tomographs. The main idea is the identification of thousands of small landmark chains that clearly correspond to the same point in a number of images. In contrast to previously published papers, our objective is to detect many short chains instead of a few long ones. This has the advantage of providing much more information to the alignment process and simplifies at the same time the task of identifying corresponding points within the tilt series. In fact, that is currently one of the limiting steps in image alignment in electron tomography: even if gold beads are available (which also result in stripe artifacts present in the 3D reconstruction that may hide interesting structural features), sometimes they are so much clustered that it is difficult to identify the same gold bead within the whole tilt series since in some of the images it is among many other gold beads clustering in the same region. Moreover, this process of identifying corresponding landmarks within the tilt series is many times done by the user, or done automatically by cross correlation of correlative image pairs and refined/corrected by the user. The main reason for the need of user interaction is that tracking the same point within the whole series is a difficult problem to be solved automatically. Placing the focus on thousands of short landmark chains makes the problem much more computationally tractable and reliable. Furthermore, the robust regression step tends to remove those possibly wrongly identified chains.
Our treatment of the chain identification is quite strict and it is one of the keys of the algorithm. The idea is to have a first rough idea of corresponding landmarks using the affine transformation between any pair of correlative images. Given the size of the images there is enough information to determine the six parameters of the affine transformation quite faithfully. Moreover, a careful global and accurate search of the parameters is performed through the combined use of a global optimizer (differential evolution) and a local optimizer (Powell's conjugate gradient method). This initial guess of the landmark chain is improved by locally shifting the image regions to maximize the correlation between corresponding regions. This allows for elastic local deformations that cannot be taken into account solely by the affine transformation (note that the elastic deformations are considered at the chain construction level, although the whole alignment procedure can only explain affine transformations and, therefore, it looses part of the richness of information captured by the landmark chains). Once the landmark chain is identified by using correlative image pairs, it is carefully refined by forcing each region to be similar not only to the equivalent regions in adjacent images, but to equivalent regions all along the image chain. This is also one of the key steps of the algorithm that greatly improve the landmark chain identification, and discards many slowly drifting chains. For the success of the algorithm it is recommended that a "dense" grid of landmark chains is available, meaning by this that there are enough landmark seeds in the image (we chose in our experiments a grid of 40 × 40) so that the local significantly overlap (in our experiments we used patches that overlapped between 35% and 70% of the patch length). We also suggest to use as small local patches as possible (values between 3% and 5% of the image size are sensible). Poor coverage of the image by the local patches will most likely result in the miss of many landmark chains, and consequently missing much of the information available in the image series.
Finally, our derivation of the equations governing the optimization of the alignment parameters, although equivalent to the previous derivation of Penczek et al. , is much more compact and provides further insight in the nature of the multivariate optimization. This fact facilitates the algorithmic implementation of the optimization procedure.
In this paper a new algorithm for automatic alignment of electron tomographs has been presented. The algorithm is based on the automatic detection of thousands of short landmark chains identifying image regions that clearly correspond to the same point in a relatively small number of images. Regions identified in this way must correlate well (usually with a correlation coefficient higher than 0.9) among all the corresponding pieces in the landmark chain, and the average of all of them. After these corresponding regions have been identified, a robust regression process is performed to finally find the relative alignment between all images. The process has been shown to successfully align images under conditions of random shifts between image pairs, and under conditions of slow but sustained drift between images. It has also been successfully applied to experimental images. The algorithm is available from the Xmipp package  and will be soon incorporated into TomoJ .
Affine registration of two images
The cross-correlation in each direction is carefully computed only in those regions which are visible in both images, i.e., only at those coordinates t H such that t H ∈ Ω i and A ij t H ∈ Ω j , being Ω i and Ω j the set of coordinates where the images i and j are defined.
The maximization problem in Eq. 2 is solved using two different optimizers: Differential Evolution (DE)  and Powell's conjugate gradient method . DE is a global optimizer based on genetic algorithms. It is rather good at approximately finding the best affine transformation between any two images. However, its convergence once it is closed to the global optimum is rather slow. Then we change to a fast local optimizer such as Powell's conjugate gradient method. In this way, we have the advantages of a global optimizer with a fast, locally convergent algorithm.
In our algorithm, image regions act as landmarks. Therefore, we must track image regions along the tilt series. However, since it is rather unlikely that the same region can be tracked along the whole series, we allow for landmark chains that start at a given image of the tilt series and finish in some other image. The region being matched along the chain need not to be detectable in all images in between the first and last image.
Every image in the tilt series, I(i), is divided in a fine grid of overlapping image regions of a size much smaller than the whole projection. The center of each region is a landmark candidate. Given a region in image i, this region is sought in image i + 1, first by predicting its position using the affine transformation Ai,i+1computed in the previous section, and then by locally optimizing the correlation index between the two corresponding regions. If the correlation index is greater than a given threshold, then the region in image i + 1 is accepted as the one matching the landmark region in image i. The new region in image i + 1 is sought in image i + 2, and the process is subsequently repeated until the correlation threshold criterium is not met. The same procedure is performed backwards. At the end of this process we have a landmark chain extending from image i - N to image i + M for some integers N and M. If the landmark chain is shorter than a certain user-supplied threshold, the chain is discarded.
If the landmark chain survives the previous filter, it is refined in an iterative way in order to avoid coordinate drifts due to the sequential search of regions. For all images between i - N + 2 and i + M the coordinate of the landmark in image l is refined with respect to that in image l - 2. This correction runs forward in the chain. When it finishes, the positions are corrected backwards (the position at image l is refined with respect to that of image l + 1) running from i + M to i - N + 1. This process is repeated with forward steps of 3, 4, ... up to a value selected by the user (i.e., the landmark at image l is refined with respect to the image l - 3, l - 4, ... respectively), while the backward correction is always performed with a step of 1. After this correction loops, the average of all regions in the chain is computed. Those regions that do not meet the correlation criterium with the average region are removed from the chain. Then the whole refinement step is repeated once again with the surviving images. Those chains whose length after refinement is smaller than the specified threshold are discarded.
Optimization of the 3D landmarks and the image in-plane rotation and shift
where V i is the set of landmarks detectable in the image i.
Derivation of the matrix form of the regression problem
A user-supplied percentage of the chains with worse fitting errors are removed from the list of landmark chains assuming that they may correspond to chains that have been wrongly computed, and all the parameters (direction of the tilt axis, position of the 3D landmarks and image in-plane rotation and shifts) are reestimated with the chains remaining. This process is repeated three times, in this way a robust fitting of the alignment parameters is performed.
On the exact location of the tilt axis
which assumes that the tilt axis passes through the coordinate system origin. However, in practice this is seldom the case, and one may wonder if it would be possible to estimate its location with respect to the microscope setup. In this appendix we show that this is not possible. This ambiguity, more than a drawback, is an advantage since it allows us centering the tilt series at will. Finally, we propose a sensible centering of the tilt axis so that the image at zero degrees (or the one closest to this tilt angle) is unmoved. However, the reader is free to locate the tilt axis to his liking.
This means that any bias we add to r0 can be automatically compensated by the appropriate shift yielding the same regression residuals. This means that absolutely locating the tilt axis within the microscope is hopeless.
This work was funded by the European Union (projects FP6-502828 and UE-512092), the 3DEM European network (LSHG-CT-2004-502828) and the ANR (PCV06-142771), the Spanish Ministerio de Educación y Ciencias (CSD2006-0023, BIO2007-67150-C01 and BIO2007-67150-C03), the Spanish Fondo de Investigación Sanitaria (04/0683) and the Comunidad de Madrid (S-GEN-0166-2006). The project described was supported by Award Number R01HL070472 from the National Heart, Lung, And Blood Institute. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Heart, Lung, And Blood Institute or the National Institutes of Health.
- Comolli LR, Kundmann M, Downing KH: Characterization of intact subcellular bodies in whole bacteria by cryo-electron tomography and spectroscopic imaging. J Microsc 2006, 223(Pt 1):40–52.View ArticlePubMedGoogle Scholar
- Donohoe BS, Mogelsvang S, Staehelin LA: Electron tomography of ER, Golgi and related membrane systems. Methods 2006, 39: 154–162.View ArticlePubMedGoogle Scholar
- Liu J, Taylor DW, Krementsova EB, Trybus KM, Taylor KA: Three dimensional structure of the myosin V inhibited state by cryoelectron tomography. Nature 2006, 442: 208–211.PubMedGoogle Scholar
- Lucic V, Forster F, Baumeister W: Structural studies by electron tomography: from cells to molecules. Annu Rev Biochem 2005, 74: 833–865.View ArticlePubMedGoogle Scholar
- Nickell S, Kofler C, Leis AP, Baumeister W: A visual approach to proteomics. Nature Reviews on Molecular Cell Biology 2006, 7: 225–230.View ArticlePubMedGoogle Scholar
- Crowther RA, Amos LA, Finch JT, De Rosier DJ, Klug A: Three dimensional reconstructions of spherical viruses by fourier synthesis from electron micrographs. Nature 1970, 226: 421–425.View ArticlePubMedGoogle Scholar
- Mastronarde DN: Dual-axis tomography: an approach with alignment methods that preserve resolution. J Struct Biol 1997, 120(3):343–352.View ArticlePubMedGoogle Scholar
- Messaoudi C, de Loubresse NG, Boudier T, Dupuis-Williams P, Marco S: Multiple-axis tomography: applications to basal bodies from Paramecium tetraurelia . Biology of the Cell 2006, 98: 415–425.View ArticlePubMedGoogle Scholar
- Penczek P, Marko M, Buttle K, Frank J: Double-tilt electron tomography. Ultramicroscopy 1995, 60: 393–410.View ArticlePubMedGoogle Scholar
- Lanzavecchia S, Cantele F, Bellon PL, Zampighi L, Kreman M, Wright E, Zampighi GA: Conical tomography of freeze-fracture replicas: a method for the study of integral membrane proteins inserted in phospholipid bilayers. J Struct Biol 2005, 149(1):87–98.View ArticlePubMedGoogle Scholar
- Zampighi GA, Zampighi L, Fain N, Wright E, Cantele F, Lanzavecchia S: Conical tomography II: A method for the study of cellular organelles in thin sections. J Struct Biol 2005, 151(3):263–274.View ArticlePubMedGoogle Scholar
- Fung JC, Liu W, de Ruitjer WJ, Chen H, Abbey CK, Sedat JW, Agard DA: Toward fully automated high-resolution electron tomography. J Struct Biol 1996, 116(1):181–189.View ArticlePubMedGoogle Scholar
- Brandt S, Heikkonen J, Engehardt P: Multiphase method for automatic alignment of transmission electron microscope images using markers. J Struct Biol 2001, 133(1):10–22.View ArticlePubMedGoogle Scholar
- Ress D, Harlow ML, Schwarz M, Marshall RM, McMahan UJ: Automatic acquisition of fiducial markers and alignment of images in tilt series for electron tomography. J Electron Microsc (Tokyo) 1999, 48(3):277–287.View ArticleGoogle Scholar
- Masich S, Östberg T, Norlén L, Shupliakov O, Daneholt B: A procedure to deposit fiducial markers on vitreous cryo-sections for cellular tomography. J Struct Biol 2006, 156(3):461–468.View ArticlePubMedGoogle Scholar
- Brandt S, Heikkonen J, Engehardt P: Automatic alignment of transmission electron microscope tilt series without fiducial markers. J Struct Biol 2001, 136(3):201–213.View ArticlePubMedGoogle Scholar
- Frank J: Computer processing of electron microscope images. Springer-Verlag 1980 chap. The role of correlation techniques in compute image processing; 187–222.Google Scholar
- Winkler H, Taylor K: Accurate marker-free alignment with simultaneous geometry determination and reconstruction of tilt series in electron tomography. Ultramicroscopy 2006, 106: 240–254.View ArticlePubMedGoogle Scholar
- Sorzano COS, Thévenaz P, Unser M: Elastic registration of biological images using vector-spline regularization. IEEE Trans Biomed Eng 2005, 52(4):652–663.View ArticlePubMedGoogle Scholar
- Scheres SHW, Núñez-Ramírez R, Sorzano COS, Carazo JM, Marabini R: Image processing for electron microscopy single-particle analysis using XMIPP. Nature Protocols 2008, 3: 977–990.PubMed CentralView ArticlePubMedGoogle Scholar
- Mastronarde DN: Electron tomography. Springer 2006 chap. Fiducial marker and hybrid alignment methods for single- and double-axis tomography; 163–186.Google Scholar
- Brandt S: Electron tomography. Springer 2006 chap. Markerless alignment in electron tomography; 187–216.Google Scholar
- Castaño-Díez D, Al-Amoudi A, Glynn AM, Seybert A, Frangakis AS: Fiducial-less alignment of cryo-sections. J Struct Biol 2007, 159(3):413–423.View ArticlePubMedGoogle Scholar
- Rieger G, Rachel R, Hermann R, Stetter KO: Ultrastructure of the hyperthermophilic archaeon Pyrodictium abyssi. J Structural Biology 1995, 115: 78–87.View ArticleGoogle Scholar
- Dubochet J, Adrian M, Chang JJ, Homo JC, Lepault J, McDowall AW, Schultz P: Cryo-electron microscopy of vitrified specimens. Quarterly Review of Biophysics 1988, 21: 129–228.View ArticleGoogle Scholar
- de Jong AF, Rees J, Busing WM, Lücken U: Integrated energy-filtering TEM – The CM120-biofilter. Prog Biophys Mol Biol 1996, 65: H116-H116.View ArticleGoogle Scholar
- Grimm R, Bärmann M, Häckl W, Typke D, Sackmann E, Baumeister W: Energy filtered electron tomography of ice-embedded actin and vesicles. Biophys J 1997, 72: 482–489.PubMed CentralView ArticlePubMedGoogle Scholar
- Dierksen K, Typke D, Hegerl R, Koster AJ, Baumeister W: Towards automatic electron tomography. Ultramicroscopy 1992, 40: 71–87.View ArticleGoogle Scholar
- Nickell S, Förster F, Linaroudis A, Del Net W, Beck F, Hegerl R, Baumeister W, Plitzko JM: TOM software toolbox: acquisition and analysis for electron tomography. J Struct Biol 2005, 149(3):227–234.View ArticlePubMedGoogle Scholar
- Sorzano COS, Marabini R, Velázquez-Muriel J, Bilbao-Castro JR, Scheres SHW, Carazo JM, Pascual-Montano A: XMIPP: A new generation of an open-source image processing package for Electron Microscopy. J Struct Biol 2004, 148(2):194–204.View ArticlePubMedGoogle Scholar
- Messaoudi C, Boudier T, Sorzano COS, Marco S: TomoJ: software for multiple-axis tomography. BMC Bioinformatics 2007, 8: 288.View ArticleGoogle Scholar
- Johnson HJ, Christensen GE: Consistent landmark and intensity-based image registration. IEEE Trans Med Imaging 2002, 21(5):450–461.View ArticlePubMedGoogle Scholar
- Arganda-Carreras I, Sorzano COS, Marabini R, Carazo JM, Ortiz-de Solórzano C, Kybic J: Consistent and elastic registration of histological sections using vector-spline regularization. Lecture Notes in Computer Science 2006, 4241: 85–95.View ArticleGoogle Scholar
- Storn R: On the use of differential evolution for function optimization. Proc NAFIPS 1996 1996, 519–523.Google Scholar
- Press W, Teukolsky SA, Vetterling WT, Flannery BP: Numerical recipes in C. 2nd edition. Cambridge University Press; 1992.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.