A Practical Approach Based on Analytic Deformable Algorithm for Scenic Image Registration

Background Image registration is to produce an entire scene by aligning all the acquired image sequences. A registration algorithm is necessary to tolerance as much as possible for intensity and geometric variation among images. However, captured image views of real scene usually produce unexpected distortions. They are generally derived from the optic characteristics of image sensors or caused by the specific scenes and objects. Methods and Findings An analytic registration algorithm considering the deformation is proposed for scenic image applications in this study. After extracting important features by the wavelet-based edge correlation method, an analytic registration approach is then proposed to achieve deformable and accurate matching of point sets. Finally, the registration accuracy is further refined to obtain subpixel precision by a feature-based Levenberg-Marquardt (FLM) method. It converges evidently faster than most other methods because of its feature-based characteristic. Conclusions We validate the performance of proposed method by testing with synthetic and real image sequences acquired by a hand-held digital still camera (DSC) and in comparison with an optical flow-based motion technique in terms of the squared sum of intensity differences (SSD) and correlation coefficient (CC). The results indicate that the proposed method is satisfactory in the registration accuracy and quality of DSC images.


Introduction
Image registration is a fundamental technology in a variety of fields and has been extensively investigated over the past few decades. It has been applied to many areas, such as medical image analysis, surveillance operations, video representation and retrieval, remote sensing, and consumer device, with different registration techniques and performance requirements [1][2][3][4][5][6][7][8]. In 2D and 3D image registration, the state-of-the-art robust techniques are reviewed and discussed regarding their advantages, drawbacks, and practical implementations [9]. Image registration in similarity measure contains feature-based and intensity-based approaches. The former accounts for important features extracted from the image. The correspondence is then established between these features by measuring the similarity. The later compare intensities or other pixel-wise signatures directly without feature extraction. The feature-based registration is effective when distinctive image features exist, while the intensity-based registration is higher computational complexity. Several well-known examples of similarity functions are the sum of squared differences (SSD), cross-correlation (CC), mutual information (MI), normalized mutual information (NMI), and Markov-Gibbs random filed (MGRF)-based. The SSD and CC are common in image registration at the same modality, while the MI and NMI are suitable in multiple modalities. However, the MI and NMI do not account for spatial relationships between adjacent pixels and they are easily influenced by image noise. MGRF-based similarity measure is derived from an MGRF model of images. The image is used as a training sample to learn a characteristic structure of pairwise pixel dependencies and Gibbs potential functions of signal co-occurrences on these pairs [9]. Image registration is mainly the process of spatially registering acquired images so that corresponding features or pixels on them are consistent in geometry. A common registration problem for the application of consumer device is to align all the acquired image sequences into a complete scene. Image alignment requires a registration algorithm that will compensate as much as possible for geometric variability among images. However, captured image views of real scene usually produce different distortions. Some are derived from the optic characteristics of image sensors, and others are caused by the specific scenes and objects. In general, we would make some reasonable assumptions to develop a fast algorithm for real time applications in consumer device fields. That is, there are no moving objects in the scenes when capturing images, and the images are acquired in short time intervals.
Another important issue for image registration is to determine the transformation model. Depending on chosen type of spatial transformation, the required number of parameters for registration model would be decided. The rigid transformation model, which preserves relative distances of points, estimates the translation and rotation, whereas the affine model [10] estimates the rigid transformation parameters and the scale factor. The affine transformation preserves collinearity. That is, parallel lines are transformed into parallel lines, and the ratios of distances are preserved along parallel lines. In addition, a more complex transformation model, perspective projection [11], takes more parameters into account. It considers not only affine transformation but also the transformations of panning and tilting. The transformation model of perspective projection is estimated to apply to the images captured from a consumer device, such as a hand-held digital still camera and a CMOS image sensor.
Moreover, a registration algorithm usually minimizes a cost function that is a combination of an objective function and smoothness constraint [12][13][14][15]. There are various algorithms that iteratively minimize the surface distance in order to linearly align two regions, such as iterative closest point algorithms.
In this study, we propose an analytic approach to achieve deformable and robust point matching. A wavelet-based method is used to extract features and discard the noise in multiscale at the same time. It then speedily evaluates the spatial correspondence and geometrical transformation between two point sets with different sizes. It is robust to noise and tolerant to distortion caused by chromatic aberration and geometry discrepancy. Finally, a feature-based Levenberg-Marquardt algorithm (FLM) is used to further refine registration results and speedily obtain subpixel accuracy because of their feature-based characteristic.
The paper is organized as follows. Section 2 addresses the proposed method in detail. In Section 3, experimental results and a discussion for some validation examples are presented. Section 4 describes the limitations of present study and future works. Finally, the conclusion is given in Section 5.

Methods
The proposed method consists of feature extraction, image registration, and registration refinement. First, feature points are extracted by wavelet-based edge correlation with large responses in multiscale. An analytic robust point matching (ARPM) algorithm is then proposed to achieve deformable and accurate registration of feature point sets. Finally, a FLM method is proposed to further obtain subpixel accuracy. The flowchart of proposed method is illustrated in Fig. 1.

Feature Extraction
It is an important issue for feature-based image registration to extract significant features from acquired images, which will produce a great influence on the registration accuracy. More specifically, feature extraction is to extract representative features from the adjacent images, so as to effectively provide the geometrical and photometric information for image registration. Multi-resolution image decomposition is a useful technique for analyzing image information at various scales. Therefore, waveletbased edge correlation, which had been verified the efficacy in feature extraction in our previous work [15], is used to extract the feature points with strong and consistent responses under different scales within a local area. In other words, feature points usually have larger values on the product of gradient moduli from multiscales, while the noise does not. In this study, Daubechies wavelet is used due to the special characteristic that it is compactly supported with extreme phase and highest number of vanishing moments for a given support width, which is beneficial to feature extraction [16].
Due to the separable characteristic of wavelet transform, we represent 2D wavelet transform as two 1D ones in x and y directions, respectively, where S(x,y)represents a 2D smoothing function. We denote w j (x,y)~1 2 2j w( x 2 j , y 2 j ) as a dilation function of w(x,y) by a scaling factor j. The gradients G j f (x,y) of an image f (x,y) in the x and y directions and its modulus M j f (x,y) at level j are described as follows, where All the edge points in image f (x,y) at level j is located with local maxima of M j f (x,y). The edge correlation, which filters out the noise by a multiscale edge confirmation to detect reliable feature points, is represented as  where n is the level number, and j is the initial level. The true feature points can be obtained by means of edge correlation. Features and noise sometimes coexist in the wavelet domain, but features can usually exist in multiscales while noise can not [15]. In addition, the level number n is chosen as two to suppress the edgebias problem of wavelets at multi-level. In this study, the observed property is used to distinguish the true feature points from noise. The procedure of feature extraction is shown in Fig. 2. A test image is given in Fig. 2(a). Fig. 2(b) shows 2D two-level wavelet decomposition for the test image. The gradient moduli of test image at level 1 and 2 are illustrated in Fig. 2(c) and 2(d), respectively. Finally, Fig. 2(e) shows the result of feature point extraction.

Image Registration
The RPM algorithm is a robust method for point-based registration [17]. Following the notation, we describe the proposed ARPM algorithm in detail.
Given two point-sets u i ,i~1,2, . . . ,H f g and v j ,j~1,2, . . . ,K È É extracted from adjacent slices, and the correspondence mapping is denoted by a matrix M consisting of m ij . The entire energy function minimized by the ARPM algorithm is as follows: where m ij [ 0,1 ½ and it subjects to P . ,H f g , T and W represent the geometric transformation and warping coefficient matrix respectively, and w stands for the warping function. The size of matrix M is Hz1 ð Þ| Kz1 ð Þ and its inner H|K portion indicates the correspondence information for two point-sets. If a point u i corresponds to a point v j , then the entry m ij of the correspondence matrix M is equal to 1; otherwise, it is assigned to zero. In addition, in order to take the outliers into account so as to still hold the constraints of the row and column summation to one, an additional row and column is appended to the suffix of the correspondence matrix M.
All the components of the energy function are interpreted in turn in the following: with the temperature parameter b is an entropy function that it avoids the elements of correspondence matrix M being negatives. The fourth term l 1 trace T{I ð Þ t T{I ð Þ Â Ã with the weighting l 1 is a constraint on the geometric transformation T by means of the penalty on the remainder of subtracting the identity matrix I from the correspondence matrix T. The final term l 2 trace W T WW ð Þ with the weighting l 2 is a constraint on the warping coefficient matrix W by means of the penalty on the warping coefficient matrix W.
As above mentioned, the minimization problem in equation (5) mainly consists of two related sub-problems: the point-sets correspondence and the geometric transformation between two adjacent slices. Given the point-sets correspondence, the geometric transformation can be evaluated by resolving the constrained leastsquares problem. Given the geometric transformation, the pointsets correspondence is found and achieved by resolving the linear assignment problem. Inspired by the idea, the algorithm incorporates the update scheme by alternating the update of the correspondence and the transformation parameters while keeping the other fixed, it is expected to jointly improve the two solutions as well as finally converge to the optimal solution.
The registration algorithm mainly consists of two principal steps. It is accomplished by using alternating update scheme. The first step is to update the point-sets correspondence matrix M as well as make sure that M corresponds to the row and column summation constraints all the time by keeping T and W fixed, with its currently evaluated transformation. Afterwards, the solution for correspondence matrix M could be calculated by means of the differentiation of the energy function in equation (5) with respect to m ij , The second step is to update the parameters of geometric transformation T and warping coefficient matrix W with the correspondence matrix M held fixed. We propose an analytic differential approach, which is namely ARPM in this study, to evaluate the parameters of T and W by means of the differentiation of the energy function in equation (5) with respect to T pq and W pq respectively, More specifically, W ij represents the element of matrix W that is located in the p th row and in the q th column. W~W ij Â Ã K|3 stands for the size of matrix W being K by 3. d½n, the unit sample sequence, is defined: d½n~1, V n~0; otherwise, d½n~0, V n=0. The two steps are iteratively performed while the temperature parameter b as well as the weighting a is gradually decreased. The decreasing process for the temperature parameter b is similar to the deterministic annealing procedure [18]. The deterministic annealing with the temperature parameter b is a procedure to adjust the flexible degree of the correspondence matrix M. The correspondence matrix M eventually approaches a binary-values matrix as the temperature b is gradually annealing. In addition, due to the property that deterministic annealing can escape from the local minima, the approach is guaranteed to obtain the near-optimal solution.

Registration Refinement
To make the registration result more accurate to further achieve subpixel precision, a refined step for image alignment is necessary. An FLM method, which is modified from [15], is proposed to apply to this study. Due to the feature-based characteristic of FLM method, it converges much faster than most other ones as long as the initial estimation is close to the global optimal solution. The sum of squared Euclidean distance errors for two feature-point sets is minimized to measure the similarity. The residue E, the sum of squared Euclidean distance errors, for the optimization criterion is written as where H and K stand for the point numbers of two point sets, which have been mentioned above. In order to resolve the nonlinear least-square minimization problem by the iterative FLM method, the Hessian matrix A and the gradient vector b must be calculated respectively, where Le ij The parameter T is then recursively updated, where l is a positive parameter, which is adjusted according to the convergence or divergence of the sum of squared errors. In other words, when the error sum decreases, the parameter l decreases and the next estimated update DT tends to the Newton method, while the parameter l increases and the next estimated update DT tends to the gradient descent approach if the error sum increases. In this study, l is initialized as 1. The FLM method is performed repeatedly until either the relative error is less than a threshold or predefined steps are reached.

Several Examples of the Registration
A hand-held digital still camera is used to capture images in the experiments. Each image is obtained with the resolution of 10246768 pixels in 24bit RGB format. A wide set of real image sequences are acquired and tested to evaluate the performance of proposed algorithm by means of the visual quality assessment. The feature extraction and registration results for a variety of image pairs are shown in Fig. 3, 4, and 5. More specifically, subfigures (a) and (b) of each figure show the image pair that will be used for registration. The results of feature point extraction from subfigures (a) and (b) are shown in subfigures (c) and (d), respectively. Finally, the registration result of image pair with proposed algorithm is shown in subfigure (e). The geometric transformation Ts of these three image pairs from Fig. 3, 4, and 5 are listed in Table 1. In addition, an optical flow-based motion algorithm [19] is implemented for comparison. The algorithm is a registration technique that takes motion estimation into account. It produces the optical flow field, a collection of two-dimensional velocity vectors, one for each small region of the image [19]. The registered results of this compared algorithm for the same image pairs are shown in Fig. 3(f)-5(f). In addition, the checkerboard visualization is also used to visualize the registration results, as shown in Fig. 3(g)-5(g), of the proposed method. The visual demonstrations indicate that the proposed method achieves better and finer results than the compared approach.
The setup of parameters in this experiment for the ARPM algorithm is described in detail as follows. The initial value for the temperature b is assigned to slightly more than the longest distance of all point pairs, and it then gradually decreases with the annealing rate 0.93. The weighting a is assigned to 5, whereas the weightings l 1 and l 2 are set to 0.1 and 0.5, respectively. The pointsets correspondence M is initialized such that all the inner entries are 1=K and the outlier ones are 1=100K . The geometric transformation T is initialized to a unit matrix. It is generally sufficient to achieve converged results by the alternating update on the correspondence M and geometric matrix T 20 runs. Fig. 3 shows the image panning and tilting problems for indoor desk and door. To show the importance of proposed registration algorithm, the image pair is captured with a large vertical motion to make big difference. Fig. 3(e) is the result of proposed registration algorithm, whereas Fig. 3(f) is the result of optical flow-based approach. As seen in these two images, it can reveal that it is difficult for optical flow-based registration approach to handle large displacement problems. The results indicate that the proposed algorithm can resolve the panning and tilting problems to achieve accurate registration, while the optical flow-based approach cannot. Fig. 4 shows the registered results of image pairs acquired from the buildings with near distances. The displacements between the image pair are quite large in this case, so the optical flow-based approach cannot obtain precise registration. The results shown in Fig. 4(e) indicate that the proposed algorithm concerns large displacements of image pairs and distortions produced from the perspective projection of acquired images. Fig. 5 shows the results of registration for image pairs with large perspective distortion in panning direction. The optical flow-based approach cannot register well since the perspective distortion is too large to accurately calculate the motion flow from matching points of image pairs. Fig. 5(e) reveals that the proposed algorithm can achieve satisfactory registration results even if the distortion in perspective projection is considerably large in titling and panning directions.

Quality of the Registration
The proposed registration algorithm has been applied to all the image pairs. To assess the quality of the registration, we calculate the mean and standard deviation of the squared sum of intensity differences (SSD) as well as the correlation coefficient (CC).

SSD~1
n where I m and I n represent a pair of images. T(.) is the geometric transform evaluated after each registration step. I m and T I n ð Þ denote the average intensities. n is the pixel number within the overlapping zone. For each image sequence, the SSD and CC provide an indirect measure of registration quality.
It is expected in theory that the image difference only shows the underlying noise from image acquisition. However, the effects of the misregistration, geometric deformation are clearly visible in registered images. In the experiments, we compare the proposed registration method with the optical flow-based motion algorithm [19]. Table 2 and 3 summarize the average results of registration quality in terms of SSD and CC for these two registration algorithms for synthetic and real DSC images, respectively. Table 2 lists the registration results for the synthetic images, which are selected from 20 DSC images with randomly selected parameters of geometric transform and Gaussian noise. The parameters for geometric transform T are selected in the ranges of , where H and W represent the height and width of images, respectively; the parameters for Gaussian noise are selected in the ranges of N(0, 5+2). Table 3 lists the registration results evaluated from all the pairs of DSC images. The results indicate that the proposed algorithm can achieve satisfactory registration accuracy and quality, which is better than the optical flow-based approach.

Statistical Evaluation
To validate whether these two algorithms are significantly different or not, one-way analysis of variance (ANOVA) is performed for the analysis of SSD and CC on both the synthetic and real data. The statistical analyses with one-way ANOVA are used to evaluate if the difference is significant for the factor SSD or CC.
We obtain p-values less than 0.0001 and less than 0.0001 for SSD and CC respectively for synthetic data, while the p-values are less than 0.0001 and less than 0.0001 for SSD and CC respectively for real data. The results of test indicate that these two algorithms are significantly different among them. More detailed comparisons  of p-values between them are then performed. The results indicate that there are significant differences in the estimation of SSD and CC between the optical flow-based approach and proposed algorithm for synthetic data (p-values be ,0.0001 and ,0.0001 for SSD and CC, respectively). In addition to synthetic data, the results of tests for real data are also discussed. The results show that the proposed algorithm is significantly better than optical flow-based approach in SSD and CC estimation (p-values be ,0.0001 and ,0.0001 for SSD and CC, respectively). Accordingly, the proposed algorithm obtains promising performance in the evaluation of registration quality for both the synthetic and real data.

Computational Cost
In addition, the computational cost is also considered to evaluate the efficiency of the proposed algorithm in this study. Image pairs are registered on a PC with an Intel Core 2 Duo E6300 processor and 1GB RAM. The results show that the proposed algorithm takes average 3.7460.25 (mean 6 standard deviation) seconds for image pairs acquired from Fig. 3, 4, and 5, whereas the optical flow-based approach takes average 5.0860.21 seconds. It indicates that the proposed method is less time consuming in computational cost.

Limitations of Present Study and Future Works
Although the proposed algorithm has less computational cost, it doesn't achieve the performance in on-line applications. In future works, the algorithm will be modified to overcome it as much as possible. In addition, the algorithm will be also applied to other style image sequences, such as medical images and satellite/remote sensing images.

Conclusion
In this study, we have presented an analytic registration algorithm for the applications of scenic images. The multi-scale concept is exploited to retain significant features and discard the noise by wavelet transform. An analytic approach is then proposed to achieve deformable and accurate registration. Finally, we refine registration accuracy to subpixel precision by the FLM method. It reduces the computational cost quite significantly due to its feature-based characteristic. It shows that this study is fairly valuable for equipping the consumers with a powerful tool in life applications.

Author Contributions
Conceived and designed the experiments: WYH. Performed the experiments: WYH. Analyzed the data: WYH. Wrote the paper: WYH. Designed the software used in analysis: WYH.