Image Alignment for Tomography Reconstruction from Synchrotron X-Ray Microscopic Images

A synchrotron X-ray microscope is a powerful imaging apparatus for taking high-resolution and high-contrast X-ray images of nanoscale objects. A sufficient number of X-ray projection images from different angles is required for constructing 3D volume images of an object. Because a synchrotron light source is immobile, a rotational object holder is required for tomography. At a resolution of 10 nm per pixel, the vibration of the holder caused by rotating the object cannot be disregarded if tomographic images are to be reconstructed accurately. This paper presents a computer method to compensate for the vibration of the rotational holder by aligning neighboring X-ray images. This alignment process involves two steps. The first step is to match the “projected feature points” in the sequence of images. The matched projected feature points in the - plane should form a set of sine-shaped loci. The second step is to fit the loci to a set of sine waves to compute the parameters required for alignment. The experimental results show that the proposed method outperforms two previously proposed methods, Xradia and SPIDER. The developed software system can be downloaded from the URL, http://www.cs.nctu.edu.tw/~chengchc/SCTA or http://goo.gl/s4AMx.


Introduction
Synchrotron X-rays possess several unique characteristics including high intensity, straight-line beam traveling, and low scattering [1], which can be used to develop an X-ray microscope [2]. Researchers have developed many synchrotron X-ray microscopy techniques for various purposes including angiography [3], X-ray lithography for nanofabrication [4], and cell tomography [5][6][7]. This study investigated a problem arising from 3D tomography reconstruction when the pixel size is in the nanoscale.
Previous research has demonstrated that a series of X-ray projections around an object can be used to reconstruct 3D volume data by using the appropriate reconstruction algorithms [8][9][10][11]. When the light source is a synchrotron X-ray, it cannot rotate around an object. Therefore, a rotatable object-holder was designed to hold an object to enable acquiring a series of projections from different angles while rotating the object (Fig. 1). The problem of mechanical imprecision arises when the resolution increases to a certain level, such as that required for cell tomography. When the pixel size is approximately 10 nm, a slight mechanical vibration can hinder accurate reconstruction. In the cell tomography experience of the authors, the pixel size of image is 11.78 nm. Rotating the object holder can cause a 5 to 30 pixel difference in position because of the mechanical instability. Although the position of the object holder can be calibrated during image acquisition, the calibration process can take an unacceptably long time, causing the object to receive excessive X-rays. The TXM controller provided by Xradia (hereinafter called ''Xradia'') was designed to solve the misalignment problem. Unfortunately, manual adjustments are generally required to obtain satisfactory tomography reconstruction. A similar problem also exists in electron microscopic tomography. Using the crosscorrelation function to align the electron microscopic images is a common solution to this problem [12]. A software system, ''SPIDER'', was implemented based on the cross-correlation function [13]. However, the cross-correlation function alignment does not considering the projection model, thus limiting the quality of tomography reconstruction.
This study presents a feature-based alignment approach for calibrating the displacement caused by mechanical vibration. Because a synchrotron X-ray is a parallel beam projection, the resulting displacement can be decomposed into vertical and horizontal displacements. The proposed method aligns the images in the vertical direction by direct image alignment. Calibrating the images in the horizontal direction is more complex than that in the vertical direction. In addition to matching feature points, the matched feature points must form sine-wave shaped loci. This study proposes approximating the loci of the matched feature points in the x-h coordinate system according to sine waves by using the least square curve fitting. The deviation between the loci and the sine waves provides information for horizontal calibration.
The remainder of this paper is organized as follows. The Results section presents the results. The Methods section presents the

Results
The proposed algorithm was used to align projection images and then reconstruct the 3D volume data of HeLa cells from X-ray projections. To verify the correctness of the alignment algorithm, a phantom data set was used to simulate the HeLa cell stained using the gold nanoparticles. The following subsections present the construction details and results obtained.

The Phantom
A volume datum containing 20 imitative gold nanoparticles was constructed, and X-ray projections of the shifted volume were generated to simulate machine vibrations. Fig. 2 shows the simulated X-ray image of the phantom cell, in which 180 projections of 5126512 pixels images were generated. The volume rotated is 1u between successive projections. Before each projection, the volume was shifted in the vertical and horizontal directions to simulate the mechanical imprecision of the object holder. The amount of displacement was determined according to a random number uniformly distributed over the range 620 pixels.
The proposed alignment algorithm was used to calibrate the images, which were first aligned in the vertical direction. The detection and matching methods of projected feature points were then applied to determine 16 feature points. The locus of the horizontal position of a feature point over various angles should form a sine-shaped curve. Fig. 3(a) shows the loci of the 16 points. This figure is in the x-h coordinate system: the vertical axis represents the rotation angle, and the horizontal axis represents the horizontal position of the projected feature point. The loci are not smooth because of horizontal displacement. The most suitable sine waves to fit the loci ( Fig. 3(b)) were then calculated, and the displacement of the feature points in the horizontal position were estimated (Fig. 3(c)). The filtered back-projection (FBP) algorithm was then used to reconstruct the 3D volume data based on the aligned X-ray images [14].
To verify the accuracy of the proposed alignment method, the positions and diameters of the spherical particles in the reconstructed images were compared with the original volume data, which included 20 spherical particles, and the diameter of each particle was 12 voxels. There were 20 particles found in the reconstructed volume. The average errors of the particle center and diameter were found to be 0.72 and 0.03 voxels, respectively. The reconstructed volume is close to the original volume. Fig. 4(a) shows a slice in the original volume data. Fig. 4(b) shows the tomography results of the same slice obtained using the proposed alignment method. For comparing the results of the proposed method with the SPIDER, the same phantom image was reconstructed using SPIDER. Fig. 4(c) shows the tomography reconstructed results using SPIDER.  The simulated X-ray image of the phantom data. 180 emulated X-ray images were generated. The size of each image is 5126512 pixels, the rotation angle between two consecutive images is 1 0 , and the ranges of the vertical and horizontal errors are +20 pixels. doi:10.1371/journal.pone.0084675.g002 The mean squared error (MSE) of the foreground voxels between the original volume data and the reconstructed volume was also calculated. A voxel was classified as a foreground voxel when its intensity exceeded a given intensity threshold. The MSE between the original volume and the reconstructed volume from the unaligned volume was found to be 0.029. The reconstructed volume obtained from the images aligned using SPIDER has an MSE of 0.002. For the volume reconstructed from the images aligned using the proposed method, the MSE was found to reduce to 0.0005. Because Xradia does not recognize the file format of the phantom images, the proposed method was not compared with Xradia.

HeLa Cells
Figs. 5(a) and (b) separately show the projection images of two HeLa cells. The first HeLa cell, named HeLa1, contained 84 identified projected feature points. Six reliable feature points were selected from the identified projected feature points. Fig. 6(a) shows the loci of the selected points in the x-h coordinate system. (c, f) were obtained using Xradia. As shown in Fig. 7, the results of the proposed method were more favorable than those of SPIDER and Xradia. Comparing the results of the SPIDER and Xradia, the proposed method exhibits most well-defined membrane structures with least amount of artifacts, which is evident from 7(a) and (d). The texture-based volume rendering algorithm [15] was used to produce a 3D image of the volume data. The proposed algorithm was applied to process 10 other X-ray image sets of HeLa cells. Among these 10 image sets, eight were successfully reconstructed (A1-A8, Figs. 9 and 10) except the other two (B1 and B2, Fig. 11). SPIDER and Xradia were also applied to the same sets of image data. Neither SPIDER nor Xradia could align the images in B1 and B2 for reconstruction. Specifically, SPIDER was effective for A3, A5, and A6; Xradia was effective for A1, A2, A5, A7, and A8. Although the reconstructions could be completed, comparing to our results, the artifacts produced due to the misalignment are more apparent. Table 1 lists the experiments conducted in this study including the phantom, HeLa1, HeLa2, and the other 10 HeLa cells, as well as the computing time required. This table shows that the main factors affecting computational time are the image size, number of images, and number of identified projected features. The experiments in this study show that when the input data contains 180 images of 102461024 pixels, the alignment can be performed in 10 min.

Sample Preparation and Image Acquisition
HeLa cells were used in this study. The cells were grown on Kapton film, and endocytosis [16] was used to stain HeLa cells by absorbing gold nanoparticles of 250 mM (micromolar). The cells were then fixed in a container using a mixture of paraformaldehyde and glutaraldehyde [17].
The synchrotron microscope used in this study was built at the National Synchrotron Radiation Research Center, Hsinchu,  Taiwan. The CCD size is 204862048 pixels, and the field of view is 24 mm. Each projection image was taken after a 1u rotation. To prevent the object holder from becoming perceptible (occurring when the object holder is nearly parallel to the X-ray), the range of the rotation angle of the object holder was 670u. Only 140 projection images were acquired. The size of each image was 102461024 pixels, and the pixel size was 11.78 nm. In this study, the exposure time of each image was 1 second.

The Alignment Method
A projected feature point is a pronounced mark in an X-ray projection image. Alignment is accomplished by first aligning the projected feature points in the vertical direction and then in the horizontal direction. The projected feature points should be maintained in the vertical direction. Thus, vertically aligning the projected feature points in the second image to the previous image is sufficient. However, the location of the feature points in the horizontal direction varies among projection images. Calculating the horizontal location of the feature points is a more difficult task than alignment in the vertical direction. The following subsections describe these steps.
Vertical Direction Alignment. For each pair of projection images, the sum of the intensity values on each row is calculated. The sums of the rows form histograms that should be similar in a pair of consecutive images. The vertical displacement can be calculated by minimizing the difference between the histograms.
Given an N6M image I(x,y), 0ƒxvN, 0ƒyvM and I(x,y)[[0, 1]. The vertical histogram h is calculated by Assume that I a is the unaligned image and I b is the reference image. The vertical correction of I a isŷ y, which can be estimated byŷ To achieve the most favorable results, the image is preprocessed to enhance the features. In this experiment, the estimated correction is more accurate when the images are enhanced by applying the edge detection method [18].
Horizontal Direction Calibration. The horizontal calibration is based on the projected feature points forming a sine-shaped locus in the x-h coordinate system. This calibration involves three   steps: detecting projected feature points, matching projected feature points to construct a set of loci from the matched projected feature points, and fitting the loci to sine curves to adjust the horizontal displacement of images.
Detecting Projected Feature Points. Feature point extraction is a fundamental step in image stitching, object recognition, and feature-based image alignment [19]. Researchers have proposed many feature detection methods. The corner detection method proposed by Harris and Stephen [20] is commonly used to extract corner-shape regions in an image. To achieve scale invariance, Kadir and Brady [21] selected the salient region from the image scale-space as the feature that possesses the maximum entropy. Lowe [22] proposed the scale-invariant feature transform (SIFT) algorithm to select local extrema from the differences of a Gaussian (DoG) pyramid in an image. The SIFT algorithm uses the gradient location-orientation histogram as a feature descriptor to achieve rotation invariance and illumination invariance. Researchers have proposed several improved versions of the SIFT algorithm. Bay et al. used the Haar wavelet to expedite feature detection [23]. Rady et al. proposed entropy-based feature detection [24], and Suri et al. combined mutual information alignment with the SIFT algorithm [25].
In this study, a modified SIFT algorithm was employed to extract automatically the projected feature points contained in Xray images. The typical SIFT implementation involves describing a feature according to its location, size, the orientation of the sampling region, and the image gradient histogram in the sampling region. Because the proposed method matches the projected feature points in two X-ray images based on mutual information [26][27][28], each projected feature point in this study contained the entropy of the sampling region rather than the image gradient histogram. To reduce the noise and the number of low-contrast projected feature points, the entropy of each selected projected feature point must exceed a given threshold. The experiments in this study entailed setting a threshold between 0.5 and 1.0. Because the features in the objects are gold nanoparticles, the size and orientation of the sampling region were fixed in this implementation.
Matching Projected Feature Points. Let F i , i~1, . . . ,m be the sets of projected feature points in m projection images. The projected feature points are classified into k groups. In the ideal case, each group is the set of projected feature points, which are the projections of a feature (i.e., gold nanoparticle) in the object from various angles. Because the rotation angle of the object is small, the projected feature points are in proximity and have similar mutual information in two consecutive images. However, the distance between the two matched projected feature points depends on the distance between the feature and the rotation axis of the object. This means that an affine transform cannot match the projected feature points in two images. Therefore, this study presents a greedy method for classifying the projected feature points. For each pair of images, the random sample consensus (RANSAC) method [29] was first applied to compute an initial alignment of the two images, and a tracking method was then employed to match the projected feature points in the next image.
Several feature tracking methods are available [19]. The proposed method is designed based on the Shafique and Shah's method [30], which is a greedy algorithm for tracking moving objects in videos, and the Tang and Tao's method [31], which integrates the hidden Markov model to eliminate unreliable matches.
Given the projected feature points sets F i{1 and F i , the RANSAC method was applied to compute a translation matrix T i,i{1 , so that a sufficient number of projected feature points p and q in F i{1 and F i respectively, jT i,i{1 (q){pj is less than a given threshold . Applying translation matrices T i,i{1 , i~2, . . . ,m to the consecutive images achieves the initial alignment of the m projection images. All of the images are aligned based on the first image.
Given the initially aligned projected feature points F i ,i~1, . . . ,m, the following procedures yield a set of possible loci of the projected feature points produced by feature points in the object. compute T i{1,i (p) where p is the final point of l and T i{1,i is the inverse of T i,i{1 . Let L be a region in I i centered at T i{1,i (p). Search in L for the projected feature points (Fig. 12). If this region contains only one projected feature point q, then that point q is selected as the final point of l. If the region contains more than one projected feature point,  3. Reverse the image orders and repeat Step 2, but do not include 2b to backtrack all loci.
The X-ray images used in this study measured 1024|1024 pixels, and 128|32 pixels was the size of search region L, and five previous points for t.
Because two loci could intersect (i.e., two projected feature points on two loci could overlap or be extremely close), the average entropy must be computed to select the best-matching projected feature point in Step 2a. In this step, some projected feature points with a high entropy in F i are not included in any locus. These significant projected feature points should not be disregarded, and Step 2b entails creating a new locus for each of them.

After
Step 2, the forward feature tracking required to construct the set of loci is complete. To verify the correctness of the loci and complete the loci added in the Step 2b, backtrack all loci in the final step.

Horizontal Displacement Estimation
Let the set of k loci collected in the previous step be fl 1 ,l 2 , . . . ,l k g. Each locus, l j , 1ƒjƒk, consists of m projected feature points, vf 1 j ,f 2 j , . . . ,f m j w. These projected feature points are expected to be the projections of a feature, p j , in the object from various angles. Assume that these projected feature points are the projections from {p=2 and p=2. Recall that f i j is not a point, but a rectangular box. Take the X -coordinates of the center of the rectangular boxes, x, and transform the j to h, the direction of the projections. Then, draw (x,h) for all f i j in the x-h coordinate system. The locus of (x,h) corresponds to a sine curve, (i.e., the sinogram). Given a set of loci of features, the horizontal alignment is conducted by fitting the curves to a set of sine curves and then by computing the deviations.
Consider a locus, l j , and assume that the projected feature points on the locus are the projections of the feature point p j . This feature point p j can be expressed as (c j ,v j ), where c is the radial coordinate and v is the angular coordinate. According to the Radon transform [8], the corresponding horizontal position x j in the i th image (rotation angle h i ), 1ƒiƒm, can be written as or in matrix form as x~Hu, where u~½u 1 ,u 2 T~½ c j sin v j ,c j cos v j T , x~½x 1 j ,x 2 j ,:::,x m j T , and  : u 1 and u 2 can be solved by the least-squares method, Finally (c j ,v j ), and The horizontal displacementx x i j is estimated bŷ Because both v j and p{v j are solutions to (7), choose the one that minimizes the sum of errors, P m i~1x x i j . To determine the horizontal displacement c i for the i th image from k feature loci, use the average of the k horizontal corrections: Some loci are unreliable because of noises, out-of-view projected feature points, or bad projected feature points matching. These loci should be removed. For each locus l j , adjust the points on l j based on the estimated c i , 1ƒiƒm. l j is unreliable if the absolute peak distance between a point on the adjusted l j and its best-fitting sine wave is greater than a given threshold s. If there are no unreliable loci, stop the algorithm and output the aligned results. Otherwise, remove the unreliable loci and repeat the horizontal displacement estimation algorithm.

Discussion
This paper presents an image alignment method for X-ray images produced by synchrotron-radiation microscopy. The proposed method enables reconstructing the 3D volume of small objects. The proposed method identifies projected feature points and classifies the projected feature points into a set of loci. The key idea of this method is that the projection of a point in the object should be a sine wave in the x-h coordinate system. Thus, fitting the set of loci to a set of sine waves can compute the parameters required for alignment.
The proposed method was applied to 12 cases of HeLa cells, and only two of the 12 could not be reconstructed. Compared with the available software systems, SPIDER and Xradia, they could respectively construct three and five cells. Neither SPIDER nor Xradia could construct the two unsuccessful cases of the proposed method. The main reason for the unsuccessful cases is the insufficient number of projected feature points in the X-ray images. The most crucial factor affecting the performance of the proposed method is the number of reliable projected feature points. If there are enough reliable projected feature points, even if the projected feature points are not in the field of view in some projections, the method still works well because it also considers the set of partial loci. The proposed method performs most favorably if the features that produce the projected feature points are close to the rotational axis. Carefully adjusting the rotation axis before images acquisition can improve the quality of the reconstruction.
Considering the shape of the projected feature points, aside from particle objects, the proposed method can manage any shape of object if it contains a sufficient number of distinct projected features. For examples, the corners of a square, the two tips of a rod, and the branch points of a tree-structured object can be used as feature points as long as the features do not deform during image acquisition. If the images satisfy these requirements, then the proposed method can successfully align the images. Table 1 lists the computing time required for the 10 successful reconstructions of HeLa cells. A graphical user interface software system for the Windows system and Mac OS X 10.8 has been developed. The software system can be downloaded from the following URL: http://www.cs.nctu.edu.tw/ cheng~chc/SCTA , or http://goo.gl/s4AMx.