Geometry-Based Global Alignment for GSMS Remote Sensing Images

Alignment of latitude and longitude for all pixels is important for geo-stationary meteorological satellite (GSMS) images. To align landmarks and non-landmarks in the GSMS images, we propose a geometry-based global alignment method. Firstly, the Global Self-consistent, Hierarchical, High-resolution Geography (GSHHG) database and GSMS images are expressed as feature maps by geometric coding. According to the geometric and gradient similarity of feature maps, initial feature matching is obtained. Then, neighborhood spatial consistency based local geometric refinement algorithm is utilized to remove outliers. Since the earth is not a standard sphere, polynomial fitting models are used to describe the global relationship between latitude, longitude and coordinates for all pixels in the GSMS images. Finally, with registered landmarks and polynomial fitting models, the latitude and longitude of each pixel in the GSMS images can be calculated. Experimental results show that the proposed method globally align the GSMS images with high accuracy, recall and significantly low computation complexity.


Introduction
In many applications, such as weather forecast, environmental monitoring and so on, determining the latitude and longitude of each pixel in the GSMS images is of great importance.However, the GSMS images have the characteristics of round-the-clock, all-weather, long range and high-resolution, which bring new challenges to practical applications.
Remote sensing images matching algorithms are usually divided into two categories: area-based methods and feature-based methods [1,2].Area-based matching algorithm establishes correspondence between two images by similarity measurements based on correlation functions.There is some classical arithmetic such as cross-correlation [3] and root mean square error (RMSE) [4].A rough-location method [5] was proposed to locate the remote image with specific physiognomy.By matching the remote sensing image and the digital map, researchers can roughly locate the remote images and the location error is less than 10 km.However, the GSMS images are generally polluted by illumination, scale variation, cloud influence and other factors, and those algorithms do not work well.A feature-based matching algorithm is widely applied to remote sensing images [6][7][8][9] because of its robustness.For example, scale-invariant feature transform (SIFT) [10,11] has an excellent performance in most circumstances.However, few feature points can be extracted from the GSMS images with SIFT due to poor textures.In addition, feature-based alignment, which only uses local gradient distribution, will lead to low precision because of too many similar features in the GSMS images.
The challenges of these points-matching methods are removing the outliers.The presence of outliers will have a negative effect on the accuracy of the matching results [12,13].To remove outliers, many algorithms based on geometric constraint and spatial information are commonly used.Among these algorithms, Random Sample Consensus (RANSAC) [14] is one of the most popular algorithms.It selects a sample randomly from the consensus set in each iteration and finds the largest consensus set to calculate the final model parameters.When the outlier is in the minority, RANSAC performs well and robustly.When the outlier is in the majority, using RANSAC will be time-consuming and unstable.By exploring the spatial relationship of matching points, a matching strategy using spatial consistent matching [8] was proposed to remove outliers.In [15][16][17], the authors proposed a spatial coding algorithm for image search, which relies on relative position relationship between pairs of matching feature points.It takes into account all matching feature pairs and encodes their coordinates to discover false matches between two images.However, the spatial relationship consistency in this method is too strict for landmark alignment.Since the earth is not a standard sphere, position deviation exists in the GSMS images.Spatial relationship consistency is effective only in a small region, and it also causes lots of correctly matched features to be deleted mistakenly.Furthermore, the number of landmarks is so large that it slows the process of removing outliers.Aguilar et al. [18] proposed a method called Graph Transformation Matching (GTM).It establishes a K-Nearest-Neighbor (KNN) graph to express neighbor geometric structures of the feature points.The mismatching feature points are determined according to the differences between KNN graph established in two images.Shi et al. [19] proposed an image registration algorithm using point structure information.After obtaining robust initial matching point pairs, the final matching results are estimated using GTM based on the local structure information of the point to remove outliers from initial correspondences.On the basis of the GTM algorithm, Weighted Graph Transformation Matching (WGTM) algorithm [20] was proposed.Utilizing the angular distances between edges that connect a feature point to its KNN as the weight, WGTM algorithms can only deal with pseudo isomorphic structures to a certain extent.This arises because angular distance is only invariant with respect to scales and rotations, and shear deformations are not considered in that case.Liu et al. [21] proposed the Restricted Spatial Order Constraints (RSOC) algorithm using a filtering strategy based on two-way geometric order constraints and two decision criteria restrictions.However, when the K-Nearest-Neighbor of the outliers are all the same, RSOC failed to remove such outliers.Zhang et al. [22] proposed a triangle-area representation of the K nearest neighbors (KNN-TAR).It utilizes the descriptor KNN-TAR to find the candidate outliers and removes the real outliers by the local structure and global information.In [23], an algorithm based on integrated spatial structure constraint (ISSC) was proposed for remote sensing image registration.First, a global structure constraint is constructed for each correspondence out of the tentative set to increase the number of inliers and raise the correct rate simultaneously.Then, a local structure constraint based on the triangle area representation is utilized on the neighboring points of each correspondence to remove outliers.Recently, Zhao et al. [24] proposed a vertex trichotomy descriptor.It utilizes the geometrical relations between any of the vertices and lines, which are constructed by mapping each vertex into trichotomy sets.A recovery and filtering vertex trichotomy matching (RFVTM) algorithm was designed to recover some inliers based on identical vertex trichotomy descriptors and restricted transformation errors.
A lot of work has been done toward the images alignment problem.Previous works can be classified in two main categories: direct [25] and feature-based methods [26,27].Direct approaches minimize pixel-to-pixel dissimilarities.While the feature-based approaches first locate a sparse set of reliable features in the image and then recover the motion parameters considering their correspondences.Miller et al. [28] proposed the congealing method by using an entropy measure to align images with respect to the distribution of the data.Cox et al. [29] proposed a least squares congealing algorithm that minimizes the sum of squared distances between images.Minimization of a log determinant cost function [30] is utilized to align images.
Inspired by these approaches, we propose a geometry-based global alignment method to align GSMS remote sensing images.According to the geometric and gradient similarity of feature maps from the GSHHG and GSMS images, initial feature matching is obtained.Then, feature refinement with a neighborhood spatial consistent matching (NSCM) algorithm is used to remove outliers.Finally, polynomial models are fitted to describe the offsets' tendency according to the matched points set.With the fitted polynomial models, the latitude and longitude of all pixels in the GSMS images can be determined.

Local Feature Matching by Geometric Coding
The shorelines of the GSHHG database correspond to the edges of the GSMS images [31], which means that shorelines can be used to simplify alignment of GSHHG and GSMS images.
Since the GSHHG database consists of polygon and line type, the size of the GSHHG database is much smaller than other reference data such as digital elevation model and digital vector map.With sub-satellite point (longitude α 0 , latitude γ 0 ) and satellite height H, the landmarks in GSHHG are mapped to a two-dimensional plane by perspective projection.Therefore, the GSHHG database is quantized to a binary image.As shown in Figure 1a, the white pixels are the landmarks defined in the GSHHG database.
The GSMS image is normalized [32,33] so that the GSHHG and GSMS images have the same size.The edges of the GSMS image extracted by Structured Forests [34] are defined as the edge probability image.As shown in Figure 1b, each element denotes the probability of the pixel being an edge candidate.To distinguish edge candidates from noise, the probability image is binarized to generate the edge binary image as depicted in Figure 1c.For a landmark in the GSHHG image, the neighborhood coding matrix W can be constructed.For a pixel P S i x S i , y S i in the edge probability image, the neighborhood coding matrix P can be constructed.Similarly the neighborhood coding matrix P can be generated with the edge binary image.The matrix W, P and P all have the same size (2K + 1) × (2K + 1).
Then, local features are matched by comparing their geometric similarity and gradient similarity.The geometric similarity between a landmark P G i x G i , y G i in the GSHHG image and a pixel P S i x S i , y S i in the edge binary image can be calculated as follows: where the W i s,t and P s,t separately denotes the s-th row and t-th column element in matrix W and P .Similarly, the gradient similarity between a landmark in the GSHHG image and a pixel P S i x S i , y S i in the edge probability image can be calculated by: The number of landmarks located within the template is calculated as follows: Both geometric and gradient similarity are measured to match local features.The procedure of local feature matching between the GSHHG and GSMS image is shown in Algorithm 1. Figure 1d shows the result of initial feature points matching.Input: W, P, P ; threshold t 1 , t 2 (t 1 is set as 0.5, t 2 is set as 0.9 based on experience) Output: the best matching pixel P S i for landmark Since there are lots of similar features in the GSMS image, local feature matching will lead to mismatching.The red circles in Figure 1d show mismatched features.The mauve circles in Figure 1d present many-to-one matched features due to the aperture effect.
The geometric relationship between matched features should not change too much across images.Based on this principle, we propose a neighborhood spatial consistent matching (NSCM) algorithm to remove outliers whose offsets between matched features have sudden mutations.
After Section 2.1, the matched set can be denoted as: M = (P G i , where the superscripts "G" and "S" refer to the GSHHG and GSMS images, respectively, (P G i , P S i ) denotes a pair of matched features and N is the number of matched features.Giving one landmark in the GSHHG image, the n nearest landmarks can be represented as The neighborhood offsets of the matched feature pair P G i , P S i can be formulated as: where ) and is constrained to ∑ µ j = 1.In addition, k is a constant normalizing µ j .When P G ij is closer to P G i , the scalar weight µ j assigns higher weights to Dx ij and Dy ij .The offsets between P G i x G i , y G i and P S i x S i , y S i in row and column can be calculated by the following formula: For the given matched feature pair P G i , P S i , the neighborhood spatial consistent matching indicates that the x i and Dx i should not deviate too much.Similarly, the y i and Dy i also should be close.This constraint can be determined: where δ and ε are two thresholds controlling sensitivity on deformations.If their values are large, the incorrect matched features are more likely to be regarded as inliers.They are both set to 0.5 according to experimental results.If P G i , P S i satisfies the low distortion constraint, it is considered as an inlier.Figure 2 is the illustration of mismatched features and many-to-one matched features.As shown in Figure 2a, (point 3, point 3') is a pair of mismatched features.The offsets between them in row and column are −2 and −2.The offsets between other pairs in neighborhood are 1 and 2. Since the offsets of (point 3, point 3') are over thresholds, they are removed.In Figure 2b, (point 2, point 2') and (point 3, point 3') are pairs of many-to-one matched features.The offsets between point 3 and point 3' in row and column are 2 and 4. The offsets between other pairs in its neighborhood are 1 and 2. (point 3, point 3') is removed and (point 2, point 2') is considered an inlier.
The details of the initial matching result and feature refinement in the southern coastal area of Thailand are shown in Figure 3. Figure 3a,c present the details of the top red circles and mauve circles, respectively, in Figure 1d.As shown in Figure 3b,d

Pixel Alignment Based on Polynomial Fitting
The earth is not a standard sphere.When using the sphere model to describe the earth, the further the pixel is away from the projection center point, the larger its distance distortion.In this case, the transformation model between sphere and plane is not suitable to describe the projection model of GSMS image.
However, the offsets between the GSHHG and GSMS images in rows and columns are smooth without distortion.For the point in the GSHHG image and its corresponding point P S i x S i , y S i in the GSMS image, the offsets between them in row and column are presented as x i and y i according to Equation (6).In order to fit the tendency of offsets in rows and columns, the polynomial functions are applied.Based on the m-th order polynomial function, the fitting functions can be defined as: where a 0 ∼ a m , c 0 ∼ c m , b 0 and d 0 are the coefficients treated as the independent variables.The point in the matched set and its corresponding x i are used to estimate the coefficients of polynomial fitting function The correlation coefficient and RMSE are considered to select the optimal coefficients.The fitted function f x i (x G i , y G i ) presents the offsets in rows changing with the coordinate (x G i , y G i ).Similarly, the coefficients of polynomial fitting function can also be estimated with the point and its corresponding y i .In addition, the fitted function f y i (x G i , y G i ) describes the offsets in columns changing with the coordinate.The offsets of pixels between the GSHHG and GSMS images can be obtained by the polynomial fitting functions in the GSHHG image, the relationship between it and its corresponding point (x S i , y S i ) in the GSMS image can be calculated as: For each pixel in the GSHHG image, the latitude and longitude information is already known.Polynomial fitting functions align all pixels of GSHHG with GSMS images globally.Therefore, the latitude and longitude of all pixels in the GSMS image can be obtained.

Dataset and Evaluation Criteria
The remote sensing images used in this experiment are from the FengyunII D meteorological satellite whose sub-satellite point is near (86 • E, 0 • N).Concerning radial distortion, only landmarks located within ±60 • of longitude and ±60 • of latitude around sub-satellite point are chosen as reference data.The size of GSMS image is normalized to 10, 000 × 10, 000 pixels.Considering efficiency, both the GSHHG and GSMS images are divided into patches [35][36][37] whose size is S1 × S2 pixels.Furthermore, feature points are matched in each pair of patches.Some shorelines can not be detected in the GSMS image due to the occlusion of clouds, causing difficulty in matching these shorelines.To reduce this difficulty, 25 patches with relatively more edges in the GSMS image are selected to perform the local feature matching and feature refinement with NSCM.
To evaluate the performance, the ground truth is manually selected from the points with the maximum gradient within their neighborhood.For each landmark in the GSHHG image, we find its corresponding point in the GSMS image as accurately as possible.Since the ground truth is labelled manually, there may be very small errors.If the distance between ground truth and matched point is no bigger than one pixel, this matched point is considered to be correct.Special attention is needed so that our manually labelled ground truth does not contain those landmarks under the clouds and fogs in the GSMS image.
In our experiments, three evaluation criteria including precision, recall and RMSE are mainly used: where N inliers represents the number of inliers in the matched set, N outliers represents the number of outliers in the matched set, N groundtruth represents the number of points of the ground truth, N p represents the number of matched pairs, P i represents the matched points and P i represents the matched points of the ground truth in the GSMS image.

Local Feature Matching by Geometric Coding
The size of the template is a key parameter for geometric coding based local feature matching.Figure 4 shows the precision and recall with K varying from 20 to 40.If the size is too small, more points are matched combined with more mismatched points.Therefore, the precision and recall are lower.As K increases, the precision is increasing and finally tends to be stable.If the size is too large, the recall is decreasing since the number of the obtained matched features is decreasing gradually.Considering the tradeoff between precision and recall, K is set to 30 in our experiments.

Feature Refinement with Neighborhood Spatial Consistent Matching (NSCM)
The NSCM algorithm is applied to remove the outliers caused by similar features and aperture effect.In the NSCM algorithm, the n nearest matched pairs are selected as neighborhood reference pairs.As depicted in Figure 5, with the value of n increasing, more neighborhood spatial consistent information is utilized and more outliers are removed.However, the spatial constraints also become stricter and the recall is decreasing.Considering the tradeoff between precision and recall, the value of n is set as 17 in the feature refinement process.

Comparison among Feature Matching Algorithms
The proposed NSCM approach is compared with seven refinement algorithms: RANSAC [14], GTM [18], WGTM [20], RSOC [21], KNN-TAR [22], ISSC [23] and RFVTM [24].Figure 6 presents the performance of these eight algorithms.In addition, the mean of experimental results are shown in Table 1.Table 1 indicates that the average precision of NSCM is the highest and the recall of NSCM algorithm ranks as medium.However, the subsequent processing can improve our recall on the basis of high precision.The RMSE value of NSCM is the smallest as shown in Table 1.As shown in Table 1, NSCM significantly outperforms the other algorithms with respect to time efficiency.Assuming that there would be N feature pairs in the matched results.In this paper, n is set to 17, which is much smaller than N. Computation complexity of NSCM is O(n × N 2 ) = O(N 2 ).

Pixel Alignment Based on Polynomial Fitting
Based on the matched set obtained by feature refinement with NSCM, the offsets between the GSHHG and GSMS images in rows and columns are fitted.The Interpolant, Lowess and Polynomial fitting types are used to get an optimal solution by comparing their precision, recall and RMSE.Table 2 shows the statistical results of the three common fitting functions.The precision of Polynomial fitting is slightly higher compared with Interpolant fitting and Lowess fitting.The recall of Polynomial fitting is far larger than the others, and the RMSE is slightly smaller than the others.In conclusion, the Polynomial fitting outperforms the other methods in all evaluation criteria.Figure 7 shows the results of Polynomial fitting functions with different order m from 1 to 5. As shown in Figure 7a,b, when m is smaller, the precision and recall are lower due to under-fitting.However, high-order polynomial leads to over-fitting.When m becomes large, the precision and recall suddenly become very low, but the RMSE becomes very high.Therefore, the third-order Polynomial fitting functions are utilized to fit the offsets' tendency.Table 3 gives the three mean values including precision, recall and RMSE before and after pixel alignment.The values of precision are close, but the recall after pixel alignment increases greatly.Figure 8 shows the result of landmark alignment.All pixels in the GSMS remote sensing image are precisely located.With pixel alignment, the latitude and longitude of all pixels in the GSMS image can be calculated.For each pixel p i , the intensity and longitude α i , latitude γ i are achieved by NSCM and Polynomial fitting.The coordinate of p i in the sub-satellite-based earth coordinate system can be represented as: where R is the radius of the earth; α 0 and γ 0 are the longitude and latitude of the sub-satellite point.With the coordinate (X i , Y i , Z i ) and intensity, the GSMS image can be displayed as a 3D earth as shown in Figure 9.

Conclusions
In this paper, we implement global alignment of all pixels in the GSMS images.Before global alignment, we do feature match between the landmarks of GSHHG and the edges of the GSMS images by geometric and gradient similarity measurement.Using spatial consistency of the matched pairs, feature refinement with a neighborhood spatial consistent matching algorithm is proposed to remove outliers.According to the experimental results, compared with other methods, our algorithm can achieve higher accuracy and lower RMSE while its time cost is significantly less than other methods.Based on polynomial fitting, global pixel alignment is applied to obtain the latitude and longitude of all pixels in the GSMS images and improve the recall significantly.The future work will focus on three-dimensional spherical stitching of multi-view remote sensing images.

Figure 1 .
Figure 1.The landmarks and GSMS images in the southern coastal area of Thailand and their initial matching results.These points in circles are outliers.(a) landmarks; (b) edge probability image; (c) edge binary image; (d) initial matching.

Figure 2 .Figure 3 .
Figure 2. Illustration of mismatched features and many-to-one matched features.(a)mismatched features; (b) many-to-one matched features.

Figure 4 .
Figure 4. Performance of local feature matching with different Ks.

Figure 5 .
Figure 5. Mean precision and recall values of feature refinement with different ns (the number of candidate matched pairs nearest the seed matched pair).(a) mean precision; and (b) mean recall.

Figure 7 .
Figure 7. Mean precision, recall and RMSE values of Polynomial fitting with different ms (the order of Polynomial function).(a) mean precision; (b) mean recall; and (c) mean RMSE.

Table 2 .
Mean precision, recall and RMSE values in Interpolant fitting, Lowess fitting and Polynomial fitting with m set to 3.

Table 3 .
Mean precision, recall and RMSE values before and after pixel alignment.