# Fast Image Registration for Spacecraft Autonomous Navigation Using Natural Landmarks

Authors: Yun-Hua Wu; Lin-Lin Ge; Feng Wang; Bing Hua; Zhi-Ming Chen; Feng Yu
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

In order to satisfy the real-time requirement of spacecraft autonomous navigation using natural landmarks, a novel algorithm called CSA-SURF (chessboard segmentation algorithm and speeded up robust features) is proposed to improve the speed without loss of repeatability performance of image registration progress. It is a combination of chessboard segmentation algorithm and SURF. Here, SURF is used to extract the features from satellite images because of its scale- and rotation-invariant properties and low computational cost. CSA is based on image segmentation technology, aiming to find representative blocks, which will be allocated to different tasks to speed up the image registration progress. To illustrate the advantages of the proposed algorithm, PCA-SURF, which is the combination of principle component analysis and SURF, is also analyzed in this paper for comparison. Furthermore, random sample consensus (RANSAC) algorithm is applied to eliminate the false matches for further accuracy improvement. The simulation results show that the proposed strategy obtains good results, especially in scaling and rotation variation. Besides, CSA-SURF decreased 50% of the time in extraction and 90% of the time in matching without losing the repeatability performance by comparing with SURF algorithm. The proposed method has been demonstrated as an alternative way for image registration of spacecraft autonomous navigation using natural landmarks.


---

## Body

## 1. Introduction

Spacecraft autonomous navigation only needs to regularly check the spacecraft working conditions, eliminating the complex navigation computing tasks, which greatly reduces the manpower and ground facility requirements and the cost of space projects [1]. Furthermore, ground stations can be destroyed during wartime. However, with the use of an autonomous navigation system, spacecraft can still work well when the ground communication is under interruption. Among many kinds of autonomous navigation systems, autonomous navigation using natural landmarks based on machine vision is a newly proposed navigation method and has potential applications for future space missions.For the landmark-based autonomous navigation, landmarks are used as reference and measuring objects. At first, the landmarks combined with position information are gathered and stored onboard the satellite. For on-orbit satellite, the camera can capture ground targets that have been stored in the satellite. During the period when the satellite runs over the target, several images from different view angles can be captured. These images are used as the inputs of image matching algorithm. Once the matching succeeds, the location information corresponding to this landmark is used to determine the position of the satellite. Then the orbit of the satellite is estimated by different positions. Figure1 shows the schematic of autonomous navigation using natural landmarks.Figure 1
Sketch of landmark-based automatic navigation.Image registration is a vital technology for landmark-based spacecraft autonomous navigation. However, autonomous navigation systems should be stable with real-time characteristics. In order to satisfy these stringent requirements, image registration progress must be fast and stable. This paper aims to propose an alternative method to improve the above performances.To improve the correctness and real-time performance in the process of image matching, Sha et al. proposed a fast matching algorithm based on imageK-gray-degree clustering, which was robust and fast under the condition of nonlinear changing of local lighting, noise, target matching of irregular shape, and even complex background [2]. However, the problem of quickly obtaining a K-degree template reflecting main features of the matching object precisely still needs to be further researched. Xu et al. proposed a novel method called DFOB for the detection, orientation computation, and description of feature points. The method was computationally efficient as it was implemented by integral images. Compared with SIFT and SURF algorithms, the computational cost of this method was much lower [3]. However, this method does not work well under large affine and perspective deformations, making it unable to perform well in wide baseline matching.Zhao et al. studied a method based on (principal component analysis) PCA [4] to speed up the image registration progress. The resulted computing time was reduced to 60% compared to single gray level normalized cross-correlation matching. He and Jiang proposed a fast image matching algorithm based on discrete Hartley transform (DHT) [5], which reduced data calculation and storage. Besides, it also improved image matching accuracy and efficiency. The PCA-based method proposed by Peng Zhao greatly improved the speed performance. However, its performance is poor when the rotation is larger than 20 deg. The image matching algorithm based on DHT also increased the image registration speed compared to traditional algorithm based on (fast Fourier transformation) FFT. But its computation time reached 18 s, which cannot satisfy the real-time requirement for natural landmark-based autonomous navigation. In order to overcome the above shortcomings, a fast image registration algorithm based on chessboard segmentation is proposed. Furthermore, RANSAC algorithm is applied to remove the error-matched key points for further repeatability improvement.The reminder of this paper is organized as follows. Firstly, Section1 introduces the purpose of this research. Then Section 2 reviews the basic theory of PCA-SURF algorithm. Later, Section 3 presents a method called chessboard segmentation algorithm to speed up image registration progress. In addition, random sample consensus algorithm is also presented in this section to obtain the registration statistical results about the number of matches and mismatches. In Section 4, the time consumptions of different methods are compared to verify the advantages of our algorithm. In Section 5, some metrics are defined to evaluate repeatability of image registration result. Besides, several tests are designed to evaluate the performance of the proposed method. Finally, Section 6 summarizes the contributions of this work.


## 2. Review of PCA-SURF

Principal component analysis (PCA) is a classical feature extraction and data representation technique widely used in the area of computer vision [6, 7]. PCA-SURF, which is a combination method of PCA and SURF [8], aims to reduce the computation by compressing the data dimension.


### 2.1. PCA-SURF Description

Suppose thatUi is a 64 × 1 description vector of a reference image I1 and Vi is a 64 × 1 description vector of an object image I2. Let X denote all SURF description vectors in the reference image and the object image:
(1)X=V1V2…VnU1U2…Um,where n and m are the number of SURF features in the reference and object images, respectively. Now, a group of orthogonal projection directions is required to be found to project X into a lower dimension space. A covariance matrix can be defined as follows:
(2)R=XXT.Then the eigenvalue and eigenvector corresponding toR are λi and μi. Actually, the value of λi is related to the amount of information in its corresponding vector μi. Larger eigenvalue means more amount of information. Then a matrix is defined
(3)C=μ1μ2…μ64.Here, it is important to satisfy thatC is an orthogonal matrix. Only with orthogonal transformation, the original vectors can be invariant to their Euclidean distance or angle with each other. What is more, in order to compress the dimension with least information loss, the following relation should be satisfied:
(4)λ1>λ2>⋯>λ64.Then with the orthogonal transformation, we have(5)Y=CTX,where Y is the feature in the new space corresponding to the original SURF feature. A new equation is defined as
(6)YT=y1y2…y64.Here, each columny1y2…y64 of Y is the principal component. Compressing the dimension of the descriptor, then y64,y63,… can be removed sequentially, because the smaller the eigenvalues are, the less the information is.


### 2.2. Application of Image Registration

In practical application, two conceptions called contribution rate and cumulative contribution rate should be defined:(7)Contributionrate=λi∑k=1nλk,Cumulativecontributionrate=∑i=1mλi∑k=1mλk.Image registration task will be completed by using the firstm principal components which occupy a large proportion of the contribution rate. In order to find the best value of the cumulative contribution rate, some simulation results are presented in the following.Figure2 shows the repeatability of PCA-SURF with different cumulative contribution rates. Different lines in the figure represent the object image with different rotations. Figure 3 presents the remaining dimensions with different cumulative contribution rates. From the above comparison, it can be seen that PCA-SURF worked poorly when the rotation is larger than 20 degrees. When the rotation is less than 20 degree, the cumulative contribution rate can be chosen as 0.95; therefore, it can largely compress the dimension of the SURF descriptor with little accuracy loss.Figure 2
Repeatability for different cumulative contribution rates.Figure 3
Remaining dimension for different cumulative contribution rates.


## 2.1. PCA-SURF Description

Suppose thatUi is a 64 × 1 description vector of a reference image I1 and Vi is a 64 × 1 description vector of an object image I2. Let X denote all SURF description vectors in the reference image and the object image:
(1)X=V1V2…VnU1U2…Um,where n and m are the number of SURF features in the reference and object images, respectively. Now, a group of orthogonal projection directions is required to be found to project X into a lower dimension space. A covariance matrix can be defined as follows:
(2)R=XXT.Then the eigenvalue and eigenvector corresponding toR are λi and μi. Actually, the value of λi is related to the amount of information in its corresponding vector μi. Larger eigenvalue means more amount of information. Then a matrix is defined
(3)C=μ1μ2…μ64.Here, it is important to satisfy thatC is an orthogonal matrix. Only with orthogonal transformation, the original vectors can be invariant to their Euclidean distance or angle with each other. What is more, in order to compress the dimension with least information loss, the following relation should be satisfied:
(4)λ1>λ2>⋯>λ64.Then with the orthogonal transformation, we have(5)Y=CTX,where Y is the feature in the new space corresponding to the original SURF feature. A new equation is defined as
(6)YT=y1y2…y64.Here, each columny1y2…y64 of Y is the principal component. Compressing the dimension of the descriptor, then y64,y63,… can be removed sequentially, because the smaller the eigenvalues are, the less the information is.


## 2.2. Application of Image Registration

In practical application, two conceptions called contribution rate and cumulative contribution rate should be defined:(7)Contributionrate=λi∑k=1nλk,Cumulativecontributionrate=∑i=1mλi∑k=1mλk.Image registration task will be completed by using the firstm principal components which occupy a large proportion of the contribution rate. In order to find the best value of the cumulative contribution rate, some simulation results are presented in the following.Figure2 shows the repeatability of PCA-SURF with different cumulative contribution rates. Different lines in the figure represent the object image with different rotations. Figure 3 presents the remaining dimensions with different cumulative contribution rates. From the above comparison, it can be seen that PCA-SURF worked poorly when the rotation is larger than 20 degrees. When the rotation is less than 20 degree, the cumulative contribution rate can be chosen as 0.95; therefore, it can largely compress the dimension of the SURF descriptor with little accuracy loss.Figure 2
Repeatability for different cumulative contribution rates.Figure 3
Remaining dimension for different cumulative contribution rates.


## 3. CSA and RANSAC

### 3.1. Chessboard Segmentation Algorithm

According to [9–11], there are three strategies that can be implemented to speed up the image registration progress. The first strategy is to decrease the dimension of the key points’ descriptors. And the second one is to decompose the task and deal them with multithreads. The third is to decrease the number of key points. Here, in this paper, the latter two methods are used to improve the speed performance.Chessboard segmentation algorithm, which combined these two ideas, can greatly improve the speed without losing performance of repeatability.Figure4 shows an island image that has been segmented into N×M parts. The number in each block represents the amount of SURF features. In order to speed up image registration progress, considering the second idea, the feature extraction tasks in each block will be allocated to different threads. Besides, by considering the third idea, only some representative blocks rather than all of them will be selected. In order to satisfy these requirements, several factors needed to be considered, for example, if the block with the most number of SURF features has been selected, then the weight of blocks around the chosen block should be decreased to make sure that regional distribution is more even. So, a block selection model should be established. The coordinate of each block is defined in Figure 4. And we suppose that Si,j is the amount of SURF features of block Di,j.Figure 4
Illustration of the number of SURF features in each block.Due to the diversity of each image, the SURF feature amounts for different images must be distinct so normalization is defined as(8)S^i,j=Si,j−SminSmax−Smin,where Smin and Smax are the minimum and maximum values, respectively, of all Si,j. And S^i,j is the normalization of Si,j. It is assumed that S^i,j is the weight of importance for each block. Here is the weight matrix W
(9)W=S^1,1…S^1,M⋮⋱⋮S^N,1…S^N,M=w1,1…w1,M⋮⋱⋮wN,1…wN,M.Local field attenuation effect simulates the phenomenon that the importance of blocks around the selected block needs to decrease. The diagram of local field attenuation effect is also demonstrated in Figure4. It is supposed that block D2,2 has been selected; then the normalized SURF feature amount S^i,j of each local field with white color should be decayed. Assuming that the attenuation threshold is T, then S^i,j is updated according to
(10)S^i,j=S^i,j−T.The entire diagram flow is presented in Figure5.Figure 5
Diagram flow of CSA.The major steps of chessboard segmentation algorithm are as follows.(1) Image Segmentation. In the first step, CSA splits the source image into N×M blocks. Simultaneously, the SURF features of each block will be extracted with the help of parallel computing technology.(2) Data Normalization. Then the amount Si,j of SURF features of each block will be normalized to establish the weight matrix W.(3) Block Selection. In each step, the block with maximum weight among the candidate set will be selected and then be removed from the candidate set.(4) Local Attenuation Effect. After the maximum weight has been selected, the weight coefficient around the selected block should minus a threshold T to simulate the local attenuation effect.Two parameters will affect the stability of CSA: the number of blocks inx and y directions. The experimental determination of the number of blocks that maximizes the stability of CSA is shown in Figure 6, which is based on an image registration task by using a collection of different natural landmarks. The terms Nx and Ny in this figure represent the number of blocks in x direction and y direction. Besides, the size of the object image that needs to be matched with the reference image is 800 × 684. The origin of an image coordinate system is located at the upper left corner of the corresponding image.Figure 6
Repeatability for combinedNx and Ny.Figure6 shows the repeatability for different combination of Nx and Ny. The red marker is the best result corresponding to 100% repeatability. According to experimental results, the CSA works better when the following equation is satisfied:
(11)NyNx≈imagewidthimageheight.Besides, CSA thresholdT represents the extent of local attenuation effect and it is also an important parameter that will affect the CSA result. Here, a solution is proposed to automatically calculate the threshold. A dataset is defined as follows:
(12)α=wi,j,i=1,…,N,j=1,…,M,where wi,j is the weight coefficient in weight matrix. Then the threshold T is equal to standard derivation of weight coefficients.
(13)T=1N×M−1∑jM∑iNwi,j−w¯2,in which w¯ is the mean of dataset α. By this equation, a suitable threshold T can be obtained.Figures7–9 present the results of CSA with different thresholds. In these figures, the selected blocks have been marked by green color. In Figure 7, the threshold T=0 means that the local attenuation effect is invalid and illustrates that these blocks whose feature number is not equal to 0 are selected. Figure 8 shows the CSA result with the threshold T=0.1 or 0.2. And Figure 9 shows the CSA result with automatic determined threshold T=0.3333. It can be seen that Figure 9 gives the most typical and minimal number of blocks with the result of automatic threshold determination.Figure 7
CSA result with thresholdT=0.Figure 8
CSA result with thresholdT=0.1 or 0.2.Figure 9
CSA result with automatic determined thresholdT=0.3333.


### 3.2. Random Sample Consensus Algorithm

Once the result of image registration is obtained, it is hard to analyze statistically the match ratio because it may have hundreds of matched key points in result. In this section, an algorithm called random sample consensus [12–14] is introduced to deal with the above problem. Besides, this algorithm also can be used to remove the wrong matched key points [15].


#### 3.2.1. Perspective Transformation

Perspective transformation is a commonly used model to represent the relationship between two images from different views.Figure10 is a demonstration of perspective transformation. P1 is a point in plane 1, and its coordinate is x,y. Q1 is the corresponding point in plane 2 with coordinate u,v. In order to get the relationship between P1 and Q1, a perspective transformation matrix H is defined as
(14)H=h1h2h3h4h5h6h7h81.Figure 10
Perspective transformation.Then the relationship betweenP1 and Q1 is
(15)λuv1=Hxy1,where λ is the scaling factor. By dividing the first and second rows with the third row, it will achieve
(16)u=h1x+h2y+h3h7x+h8y+1,v=h4x+h5y+h6h7x+h8y+1.Then we assume thatH^=h1…h8T; to transfer the above equation into the matrix form, we can obtain
(17)xy1000−xu−yu000xy1−xv−yvH^=uv.In order to solveH^, four or more corresponding pairs are required. So, we suppose that there are four known corresponding groups P1,P2,P3,P4 and Q1,Q2,Q3,Q4, and satisfying
(18)x1y11000−x1u1−y1u1000x1y11−x1v1−y1v1…x4y41000−x4u4−y4u4000x4y41−x4v4−y4v4h1h2⋮h7h8=u1v1⋮u4v4,it can be simplified as
(19)AH^=b.Because of distortion, noise, or other reasons, the above equation may be a contradictory equation. Therefore, the least squares method is applied to find the satisfied result(20)H^=ATA−1ATb.Eventually, the perspective transformation matrixH can be recovered from H^.Assuming thatPx,y is the key point in the first picture. Qu,v is the corresponding matched key point in the second image. H1 is the perspective transformation matrix between the first and second image. Then we have
(21)λx′y′1T=H1xy1T,where x′,y′ is the corresponding point after perspective projection transformation on key point P. Then the perspective projection error can be defined as follows:
(22)Err=x′−u2+y′−v2.


#### 3.2.2. RANSAC Algorithm Progress

Figure11 shows the diagram flow of RANSAC algorithm. The main stages of RANSAC algorithm are described as follows:(1) Initialization. The threshold of minimum perspective project error Tppe should be initialized at first. Besides, the maximum iterator time Miter should also be initialized here.(2) Perspective Transformation Solving. The perspective transformation H can be solved by least squares method with four groups of key points, among which four key points are randomly selected from the dataset of source image key points. The other four key points are randomly selected from the dataset of target image key points.(3) Perspective Matrix Validation. Once a new perspective projection error has been obtained in each iterator, a validation needs to be performed. The minimum of perspective projection error Err between the key points of the source image and the key points of the target image will be calculated. If Err is smaller than the smallest perspective projection error Emin which is stored before, the minimum perspective projection error Emin will be updated.(4) Loop Termination. If Emin>Tppe and iterator>Miter, the loop will go back to step 2; otherwise, the loop will be terminated.Figure 11
Diagram flow of RANSAC algorithm.


## 3.1. Chessboard Segmentation Algorithm

According to [9–11], there are three strategies that can be implemented to speed up the image registration progress. The first strategy is to decrease the dimension of the key points’ descriptors. And the second one is to decompose the task and deal them with multithreads. The third is to decrease the number of key points. Here, in this paper, the latter two methods are used to improve the speed performance.Chessboard segmentation algorithm, which combined these two ideas, can greatly improve the speed without losing performance of repeatability.Figure4 shows an island image that has been segmented into N×M parts. The number in each block represents the amount of SURF features. In order to speed up image registration progress, considering the second idea, the feature extraction tasks in each block will be allocated to different threads. Besides, by considering the third idea, only some representative blocks rather than all of them will be selected. In order to satisfy these requirements, several factors needed to be considered, for example, if the block with the most number of SURF features has been selected, then the weight of blocks around the chosen block should be decreased to make sure that regional distribution is more even. So, a block selection model should be established. The coordinate of each block is defined in Figure 4. And we suppose that Si,j is the amount of SURF features of block Di,j.Figure 4
Illustration of the number of SURF features in each block.Due to the diversity of each image, the SURF feature amounts for different images must be distinct so normalization is defined as(8)S^i,j=Si,j−SminSmax−Smin,where Smin and Smax are the minimum and maximum values, respectively, of all Si,j. And S^i,j is the normalization of Si,j. It is assumed that S^i,j is the weight of importance for each block. Here is the weight matrix W
(9)W=S^1,1…S^1,M⋮⋱⋮S^N,1…S^N,M=w1,1…w1,M⋮⋱⋮wN,1…wN,M.Local field attenuation effect simulates the phenomenon that the importance of blocks around the selected block needs to decrease. The diagram of local field attenuation effect is also demonstrated in Figure4. It is supposed that block D2,2 has been selected; then the normalized SURF feature amount S^i,j of each local field with white color should be decayed. Assuming that the attenuation threshold is T, then S^i,j is updated according to
(10)S^i,j=S^i,j−T.The entire diagram flow is presented in Figure5.Figure 5
Diagram flow of CSA.The major steps of chessboard segmentation algorithm are as follows.(1) Image Segmentation. In the first step, CSA splits the source image into N×M blocks. Simultaneously, the SURF features of each block will be extracted with the help of parallel computing technology.(2) Data Normalization. Then the amount Si,j of SURF features of each block will be normalized to establish the weight matrix W.(3) Block Selection. In each step, the block with maximum weight among the candidate set will be selected and then be removed from the candidate set.(4) Local Attenuation Effect. After the maximum weight has been selected, the weight coefficient around the selected block should minus a threshold T to simulate the local attenuation effect.Two parameters will affect the stability of CSA: the number of blocks inx and y directions. The experimental determination of the number of blocks that maximizes the stability of CSA is shown in Figure 6, which is based on an image registration task by using a collection of different natural landmarks. The terms Nx and Ny in this figure represent the number of blocks in x direction and y direction. Besides, the size of the object image that needs to be matched with the reference image is 800 × 684. The origin of an image coordinate system is located at the upper left corner of the corresponding image.Figure 6
Repeatability for combinedNx and Ny.Figure6 shows the repeatability for different combination of Nx and Ny. The red marker is the best result corresponding to 100% repeatability. According to experimental results, the CSA works better when the following equation is satisfied:
(11)NyNx≈imagewidthimageheight.Besides, CSA thresholdT represents the extent of local attenuation effect and it is also an important parameter that will affect the CSA result. Here, a solution is proposed to automatically calculate the threshold. A dataset is defined as follows:
(12)α=wi,j,i=1,…,N,j=1,…,M,where wi,j is the weight coefficient in weight matrix. Then the threshold T is equal to standard derivation of weight coefficients.
(13)T=1N×M−1∑jM∑iNwi,j−w¯2,in which w¯ is the mean of dataset α. By this equation, a suitable threshold T can be obtained.Figures7–9 present the results of CSA with different thresholds. In these figures, the selected blocks have been marked by green color. In Figure 7, the threshold T=0 means that the local attenuation effect is invalid and illustrates that these blocks whose feature number is not equal to 0 are selected. Figure 8 shows the CSA result with the threshold T=0.1 or 0.2. And Figure 9 shows the CSA result with automatic determined threshold T=0.3333. It can be seen that Figure 9 gives the most typical and minimal number of blocks with the result of automatic threshold determination.Figure 7
CSA result with thresholdT=0.Figure 8
CSA result with thresholdT=0.1 or 0.2.Figure 9
CSA result with automatic determined thresholdT=0.3333.


## 3.2. Random Sample Consensus Algorithm

Once the result of image registration is obtained, it is hard to analyze statistically the match ratio because it may have hundreds of matched key points in result. In this section, an algorithm called random sample consensus [12–14] is introduced to deal with the above problem. Besides, this algorithm also can be used to remove the wrong matched key points [15].


### 3.2.1. Perspective Transformation

Perspective transformation is a commonly used model to represent the relationship between two images from different views.Figure10 is a demonstration of perspective transformation. P1 is a point in plane 1, and its coordinate is x,y. Q1 is the corresponding point in plane 2 with coordinate u,v. In order to get the relationship between P1 and Q1, a perspective transformation matrix H is defined as
(14)H=h1h2h3h4h5h6h7h81.Figure 10
Perspective transformation.Then the relationship betweenP1 and Q1 is
(15)λuv1=Hxy1,where λ is the scaling factor. By dividing the first and second rows with the third row, it will achieve
(16)u=h1x+h2y+h3h7x+h8y+1,v=h4x+h5y+h6h7x+h8y+1.Then we assume thatH^=h1…h8T; to transfer the above equation into the matrix form, we can obtain
(17)xy1000−xu−yu000xy1−xv−yvH^=uv.In order to solveH^, four or more corresponding pairs are required. So, we suppose that there are four known corresponding groups P1,P2,P3,P4 and Q1,Q2,Q3,Q4, and satisfying
(18)x1y11000−x1u1−y1u1000x1y11−x1v1−y1v1…x4y41000−x4u4−y4u4000x4y41−x4v4−y4v4h1h2⋮h7h8=u1v1⋮u4v4,it can be simplified as
(19)AH^=b.Because of distortion, noise, or other reasons, the above equation may be a contradictory equation. Therefore, the least squares method is applied to find the satisfied result(20)H^=ATA−1ATb.Eventually, the perspective transformation matrixH can be recovered from H^.Assuming thatPx,y is the key point in the first picture. Qu,v is the corresponding matched key point in the second image. H1 is the perspective transformation matrix between the first and second image. Then we have
(21)λx′y′1T=H1xy1T,where x′,y′ is the corresponding point after perspective projection transformation on key point P. Then the perspective projection error can be defined as follows:
(22)Err=x′−u2+y′−v2.


### 3.2.2. RANSAC Algorithm Progress

Figure11 shows the diagram flow of RANSAC algorithm. The main stages of RANSAC algorithm are described as follows:(1) Initialization. The threshold of minimum perspective project error Tppe should be initialized at first. Besides, the maximum iterator time Miter should also be initialized here.(2) Perspective Transformation Solving. The perspective transformation H can be solved by least squares method with four groups of key points, among which four key points are randomly selected from the dataset of source image key points. The other four key points are randomly selected from the dataset of target image key points.(3) Perspective Matrix Validation. Once a new perspective projection error has been obtained in each iterator, a validation needs to be performed. The minimum of perspective projection error Err between the key points of the source image and the key points of the target image will be calculated. If Err is smaller than the smallest perspective projection error Emin which is stored before, the minimum perspective projection error Emin will be updated.(4) Loop Termination. If Emin>Tppe and iterator>Miter, the loop will go back to step 2; otherwise, the loop will be terminated.Figure 11
Diagram flow of RANSAC algorithm.


## 3.2.1. Perspective Transformation

Perspective transformation is a commonly used model to represent the relationship between two images from different views.Figure10 is a demonstration of perspective transformation. P1 is a point in plane 1, and its coordinate is x,y. Q1 is the corresponding point in plane 2 with coordinate u,v. In order to get the relationship between P1 and Q1, a perspective transformation matrix H is defined as
(14)H=h1h2h3h4h5h6h7h81.Figure 10
Perspective transformation.Then the relationship betweenP1 and Q1 is
(15)λuv1=Hxy1,where λ is the scaling factor. By dividing the first and second rows with the third row, it will achieve
(16)u=h1x+h2y+h3h7x+h8y+1,v=h4x+h5y+h6h7x+h8y+1.Then we assume thatH^=h1…h8T; to transfer the above equation into the matrix form, we can obtain
(17)xy1000−xu−yu000xy1−xv−yvH^=uv.In order to solveH^, four or more corresponding pairs are required. So, we suppose that there are four known corresponding groups P1,P2,P3,P4 and Q1,Q2,Q3,Q4, and satisfying
(18)x1y11000−x1u1−y1u1000x1y11−x1v1−y1v1…x4y41000−x4u4−y4u4000x4y41−x4v4−y4v4h1h2⋮h7h8=u1v1⋮u4v4,it can be simplified as
(19)AH^=b.Because of distortion, noise, or other reasons, the above equation may be a contradictory equation. Therefore, the least squares method is applied to find the satisfied result(20)H^=ATA−1ATb.Eventually, the perspective transformation matrixH can be recovered from H^.Assuming thatPx,y is the key point in the first picture. Qu,v is the corresponding matched key point in the second image. H1 is the perspective transformation matrix between the first and second image. Then we have
(21)λx′y′1T=H1xy1T,where x′,y′ is the corresponding point after perspective projection transformation on key point P. Then the perspective projection error can be defined as follows:
(22)Err=x′−u2+y′−v2.


## 3.2.2. RANSAC Algorithm Progress

Figure11 shows the diagram flow of RANSAC algorithm. The main stages of RANSAC algorithm are described as follows:(1) Initialization. The threshold of minimum perspective project error Tppe should be initialized at first. Besides, the maximum iterator time Miter should also be initialized here.(2) Perspective Transformation Solving. The perspective transformation H can be solved by least squares method with four groups of key points, among which four key points are randomly selected from the dataset of source image key points. The other four key points are randomly selected from the dataset of target image key points.(3) Perspective Matrix Validation. Once a new perspective projection error has been obtained in each iterator, a validation needs to be performed. The minimum of perspective projection error Err between the key points of the source image and the key points of the target image will be calculated. If Err is smaller than the smallest perspective projection error Emin which is stored before, the minimum perspective projection error Emin will be updated.(4) Loop Termination. If Emin>Tppe and iterator>Miter, the loop will go back to step 2; otherwise, the loop will be terminated.Figure 11
Diagram flow of RANSAC algorithm.


## 4. Time Consumption Comparison

In order to verify the advantage of the proposed CSA-SURF algorithm, the SURF and PCA-SURF algorithms are used for comparison. Feature extraction and matching progress are different for the above three algorithms, of which the time complexities are discussed in the following.


### 4.1. Time Consumption for Feature Extraction

It is well known that SIFT or SURF algorithms use the sliding window method to detect local extrema. PCA-SURF and CSA-SURF, having the same mechanism, are the improved algorithms based on SURF. However, there are some differences among them. In order to descript more clearly, the width and height of the image are defined asIw and Ih and the step size is σ. The time for solving the eigenvalue and eigenvector is Teig. The time for transferring the old feature space to the new feature space is Ttran. The time for computing the local feature is Tlocal. The extraction times of SURF, PCA-SURF, and CSA-SURF are Tsurfe, Tpca−surfe, and Tcsa−surfe, respectively. The pseudocodes for the feature extraction progress of SURF, PCA-SURF, and CSA-SURF algorithms are presented as in Pseudocodes 1, 2, and 3.Pseudocode 1: Pseudocode for SURF feature extraction progress.
fori=1:σ:Iw−3
forj=1:σ:Ih−3
compute local feature for each pixel
end
endPseudocode 2: Pseudocode for PCA-SURF feature extraction progress.
fori=1:σ:Iw−3
forj=1:σ:Ih−3
compute local feature for each pixel
end
end
solve for eigenvalue and eigenvector
transfer feature spacePseudocode 3: Pseudocode for CSA-SURF feature extraction progress.
fori=1:n
forj=1:m
fork1=1:Iw/n−3
fork2=Ih/m−3
compute local feature for each pixel
end
end
end
endThe extraction time of SURF, PCA-SURF, and CSA-SURF can be calculated by(23)Tsurfe=Iw−3σIh−3σTlocal,Tpca−surfe=Iw−3σIh−3σTlocal+Teig+Ttran,Tcsa−surfe=nmIw/n−3σIh/m−3σTlocal,satisfying
(24)Tcsa−surfe<Tsurfe<Tpca−surfe,in which the symbol ⋅ is a rounding down function.The relationship of (24) can be proven as follows.Obviously, we have(25)nmIw/n−3σIh/m−3σ<nmIw/n−3Ih/m−3σ2,(26)nmIw/n−3Ih/m−3σ2=Iw−3nIh−3mσ2,(27)Iw−3nIh−3mσ2<Iw−3σIh−3σ.Combining (26) and (27), we obtain
(28)nmIw/n−3σIh/m−3σ<Iw−3σIh−3σ,that is,(29)Tcsa−surfe<Tsurfe.Obviously,(30)Tsurfe<Tpca−surfe.Therefore, (24) has been verified.


### 4.2. Time Consumption of Matching

In CSA-SURF, the image is divided inton blocks in horizontal and m blocks in vertical. The number of extracted features for SURF and PCA-SURF is N, while the corresponding number for CSA-SURF is M and is much smaller than N. Besides, the dimension of PCA-SURF is l. The time consumption for two features matching is Tmatch. The matching time for SURF, PCA-SURF, and CSA-SURF are Tsurfm, Tpca−surfm, and Tcsa−surfm and can be calculated by
(31)Tsurfm=N2Tmatch,Tpca−surfm=N2l64Tmatch,Tcsa−surfm=M2Tmatch.AsM is much less than N, then we have
(32)Tcsa−surfm<Tpca−surfm,whenM≪N,which can be proven as follows.As we can see,(33)M2<N2,N2l64<N2.Therefore, we have(34)Tpca−surfm<Tsurfm,Tcsa−surfm<Tsurfm.IfM≪N, then
(35)MN2⟶smallnumber.For the spacecraft autonomous navigation problem using natural landmarks, only 5 to 10 features satisfying certain relative distance constraints are required; therefore, it is easy to guarantee thatM/N2<N2/l and then we have Tcsa−surfm<Tpca−surfm.From the above comparison, the proposed CSA-SURF method consumes the least time.


## 4.1. Time Consumption for Feature Extraction

It is well known that SIFT or SURF algorithms use the sliding window method to detect local extrema. PCA-SURF and CSA-SURF, having the same mechanism, are the improved algorithms based on SURF. However, there are some differences among them. In order to descript more clearly, the width and height of the image are defined asIw and Ih and the step size is σ. The time for solving the eigenvalue and eigenvector is Teig. The time for transferring the old feature space to the new feature space is Ttran. The time for computing the local feature is Tlocal. The extraction times of SURF, PCA-SURF, and CSA-SURF are Tsurfe, Tpca−surfe, and Tcsa−surfe, respectively. The pseudocodes for the feature extraction progress of SURF, PCA-SURF, and CSA-SURF algorithms are presented as in Pseudocodes 1, 2, and 3.Pseudocode 1: Pseudocode for SURF feature extraction progress.
fori=1:σ:Iw−3
forj=1:σ:Ih−3
compute local feature for each pixel
end
endPseudocode 2: Pseudocode for PCA-SURF feature extraction progress.
fori=1:σ:Iw−3
forj=1:σ:Ih−3
compute local feature for each pixel
end
end
solve for eigenvalue and eigenvector
transfer feature spacePseudocode 3: Pseudocode for CSA-SURF feature extraction progress.
fori=1:n
forj=1:m
fork1=1:Iw/n−3
fork2=Ih/m−3
compute local feature for each pixel
end
end
end
endThe extraction time of SURF, PCA-SURF, and CSA-SURF can be calculated by(23)Tsurfe=Iw−3σIh−3σTlocal,Tpca−surfe=Iw−3σIh−3σTlocal+Teig+Ttran,Tcsa−surfe=nmIw/n−3σIh/m−3σTlocal,satisfying
(24)Tcsa−surfe<Tsurfe<Tpca−surfe,in which the symbol ⋅ is a rounding down function.The relationship of (24) can be proven as follows.Obviously, we have(25)nmIw/n−3σIh/m−3σ<nmIw/n−3Ih/m−3σ2,(26)nmIw/n−3Ih/m−3σ2=Iw−3nIh−3mσ2,(27)Iw−3nIh−3mσ2<Iw−3σIh−3σ.Combining (26) and (27), we obtain
(28)nmIw/n−3σIh/m−3σ<Iw−3σIh−3σ,that is,(29)Tcsa−surfe<Tsurfe.Obviously,(30)Tsurfe<Tpca−surfe.Therefore, (24) has been verified.


## 4.2. Time Consumption of Matching

In CSA-SURF, the image is divided inton blocks in horizontal and m blocks in vertical. The number of extracted features for SURF and PCA-SURF is N, while the corresponding number for CSA-SURF is M and is much smaller than N. Besides, the dimension of PCA-SURF is l. The time consumption for two features matching is Tmatch. The matching time for SURF, PCA-SURF, and CSA-SURF are Tsurfm, Tpca−surfm, and Tcsa−surfm and can be calculated by
(31)Tsurfm=N2Tmatch,Tpca−surfm=N2l64Tmatch,Tcsa−surfm=M2Tmatch.AsM is much less than N, then we have
(32)Tcsa−surfm<Tpca−surfm,whenM≪N,which can be proven as follows.As we can see,(33)M2<N2,N2l64<N2.Therefore, we have(34)Tpca−surfm<Tsurfm,Tcsa−surfm<Tsurfm.IfM≪N, then
(35)MN2⟶smallnumber.For the spacecraft autonomous navigation problem using natural landmarks, only 5 to 10 features satisfying certain relative distance constraints are required; therefore, it is easy to guarantee thatM/N2<N2/l and then we have Tcsa−surfm<Tpca−surfm.From the above comparison, the proposed CSA-SURF method consumes the least time.


## 5. Evaluation

### 5.1. Evaluation Metrics

In order to descript the repeatability performance, the following evaluation metrics are defined in this section. Positive in this paper means the key points that can be matched. Negative means the key points that cannot be matched.Then we suppose thatγTP is the number of positive key points that are correct matches. γFP is the number of negative key points that are false matches. γTN is the number of negative key points that correct mismatches. And γFN is the number of negative key points that are false mismatches. To explain it more clearly, Figure 12 shows the relationship of different symbols.Figure 12
Relationship of different symbols.Then we have the following definitions:(36)Precision=γTPγTP+γFP=correctmatchespositive,(37)Recall=γTPγTP+γFN=correctmatchespositive,Precision means the proportion of correct matches to positive. Recall means the proportion of correct matches to matches.Recall and precision are conflicting, so a recall versus precision graph will help in analyzing the repeatability performance.


### 5.2. Speed Test

The following experiments are tested on Core I5 2.3 GHz CPU, 6GB RAM laptop with Windows 7 operation system. VS2013 and OpenCV2.4 are used to carry out all the experiments. And all datasets are from an open source database called UC Merced Land Use Dataset. And the resolution is 256 × 256.Figures13–15 are the matching comparison results of fifty images about SURF, PCA-SURF, and CSA-SURF. It can be seen that PCA-SURF consumes much more time than SURF and CSA-SURF during feature extraction process. The CSA-SURF algorithm works faster than SURF and PCA-SURF algorithms for both feature extraction and matching processes.Figure 13
Time consumption for SURF.Figure 14
Time consumption for PCA-SURF.Figure 15
Time consumption for CSA-SURF.


### 5.3. Repeatability Test

Figure16 shows the repeatability of different methods for different conditions. But recall and precision may be conflict in some condition, so a recall versus precision graph has been established to help us analyze the repeatability performance. The red line is the test under scaling; it can be seen that CSA-SURF performs worse when precision<0.2 and its performance is not stable when 0.2<precision<0.45. For precision<0.45, it works well with recall=1. The green line represents the test of CSA-SURF algorithm under rotation. Its performance is much the same as the solid line’s case, while its unstable region is 0.2<precision<0.4. The blue line represents the test of CSA-SURF algorithm under salt and pepper noise [16] with a noise density equal to 0.35. It can be seen that CSA-SURF algorithm is not stable under salt and pepper noise, but recall still maintains its value above 0.7.Figure 16
Recall versus precision graph of three methods under different conditions.Furthermore, the above dataset with fifty images is also used to test the repeatability performance of the three methods. Figure17 is the recall of different methods for these images. The results show that the repeatability of SURF and CSA-SURF is very close. Therefore, CSA-SURF can speed up the registration progress without losing its repeatability performance. However, PCA-SURF is unstable as it can get the correct matches in some conditions.Figure 17
Recall of different methods.


### 5.4. Application on Natural Landmark Registration

The major purpose of this paper is to apply the proposed method into natural landmark registration task. Five different kinds of natural landmarks, including island, airport, railway, river, and coastline, are collected for the demonstration testing.Figure18 presents the registration results for Hawaiian islands with different algorithms, in which SURF and CSA-SURF work fine but PCA-SURF has a wrong match. Figure 19 shows the matching results of two O’Hare International Airport images with different algorithms, indicating that the PCA-SURF algorithm is unstable as it cannot find any matched key points in this example. Figures 20–22 are the matching results of Qingzang railway, river, and coastline images with different algorithms.Figure 18
Hawaiian island registration result with SURF, PCA-SURF, and CSA-SURF.Figure 19
O’Hare International Airport registration result with SURF, PCA-SURF, and CSA-SURF.Figure 20
Qingzang railway registration result with SURF, PCA-SURF, and CSA-SURF.Figure 21
River landmark registration result with SURF, PCA-SURF, and CSA-SURF.Figure 22
Coastline landmark registration result with SURF, PCA-SURF, and CSA-SURF.It can be seen from the above testing results that CSA-SURF algorithm works well for different images with significant contour features. The proposed CSA-SURF algorithm, which is improved from SURF, keeps the repeatability of SURF unchanged and improves the image matching speed. Many other pictures have also been tested; however, the results are not presented for the page limitation of the paper.


## 5.1. Evaluation Metrics

In order to descript the repeatability performance, the following evaluation metrics are defined in this section. Positive in this paper means the key points that can be matched. Negative means the key points that cannot be matched.Then we suppose thatγTP is the number of positive key points that are correct matches. γFP is the number of negative key points that are false matches. γTN is the number of negative key points that correct mismatches. And γFN is the number of negative key points that are false mismatches. To explain it more clearly, Figure 12 shows the relationship of different symbols.Figure 12
Relationship of different symbols.Then we have the following definitions:(36)Precision=γTPγTP+γFP=correctmatchespositive,(37)Recall=γTPγTP+γFN=correctmatchespositive,Precision means the proportion of correct matches to positive. Recall means the proportion of correct matches to matches.Recall and precision are conflicting, so a recall versus precision graph will help in analyzing the repeatability performance.


## 5.2. Speed Test

The following experiments are tested on Core I5 2.3 GHz CPU, 6GB RAM laptop with Windows 7 operation system. VS2013 and OpenCV2.4 are used to carry out all the experiments. And all datasets are from an open source database called UC Merced Land Use Dataset. And the resolution is 256 × 256.Figures13–15 are the matching comparison results of fifty images about SURF, PCA-SURF, and CSA-SURF. It can be seen that PCA-SURF consumes much more time than SURF and CSA-SURF during feature extraction process. The CSA-SURF algorithm works faster than SURF and PCA-SURF algorithms for both feature extraction and matching processes.Figure 13
Time consumption for SURF.Figure 14
Time consumption for PCA-SURF.Figure 15
Time consumption for CSA-SURF.


## 5.3. Repeatability Test

Figure16 shows the repeatability of different methods for different conditions. But recall and precision may be conflict in some condition, so a recall versus precision graph has been established to help us analyze the repeatability performance. The red line is the test under scaling; it can be seen that CSA-SURF performs worse when precision<0.2 and its performance is not stable when 0.2<precision<0.45. For precision<0.45, it works well with recall=1. The green line represents the test of CSA-SURF algorithm under rotation. Its performance is much the same as the solid line’s case, while its unstable region is 0.2<precision<0.4. The blue line represents the test of CSA-SURF algorithm under salt and pepper noise [16] with a noise density equal to 0.35. It can be seen that CSA-SURF algorithm is not stable under salt and pepper noise, but recall still maintains its value above 0.7.Figure 16
Recall versus precision graph of three methods under different conditions.Furthermore, the above dataset with fifty images is also used to test the repeatability performance of the three methods. Figure17 is the recall of different methods for these images. The results show that the repeatability of SURF and CSA-SURF is very close. Therefore, CSA-SURF can speed up the registration progress without losing its repeatability performance. However, PCA-SURF is unstable as it can get the correct matches in some conditions.Figure 17
Recall of different methods.


## 5.4. Application on Natural Landmark Registration

The major purpose of this paper is to apply the proposed method into natural landmark registration task. Five different kinds of natural landmarks, including island, airport, railway, river, and coastline, are collected for the demonstration testing.Figure18 presents the registration results for Hawaiian islands with different algorithms, in which SURF and CSA-SURF work fine but PCA-SURF has a wrong match. Figure 19 shows the matching results of two O’Hare International Airport images with different algorithms, indicating that the PCA-SURF algorithm is unstable as it cannot find any matched key points in this example. Figures 20–22 are the matching results of Qingzang railway, river, and coastline images with different algorithms.Figure 18
Hawaiian island registration result with SURF, PCA-SURF, and CSA-SURF.Figure 19
O’Hare International Airport registration result with SURF, PCA-SURF, and CSA-SURF.Figure 20
Qingzang railway registration result with SURF, PCA-SURF, and CSA-SURF.Figure 21
River landmark registration result with SURF, PCA-SURF, and CSA-SURF.Figure 22
Coastline landmark registration result with SURF, PCA-SURF, and CSA-SURF.It can be seen from the above testing results that CSA-SURF algorithm works well for different images with significant contour features. The proposed CSA-SURF algorithm, which is improved from SURF, keeps the repeatability of SURF unchanged and improves the image matching speed. Many other pictures have also been tested; however, the results are not presented for the page limitation of the paper.


## 6. Conclusion

Because of the restricted computational resource and horrible environment for an autonomous navigation system, traditional image registration methods cannot satisfy practical requirements in speed. A novel algorithm called chessboard segmentation algorithm is proposed to solve the above problem. Because the new method is based on SURF features, it inherits lots of their advantages, such as scale and rotation invariant properties. To verify the improvement of the proposed algorithm, the PCA-SURF algorithm which was proposed in recent years is also presented in this paper for comparison. Besides, RANSAC algorithm is applied to remove the false negative key points to further improve the accuracy of the proposed algorithm. Thorough experiments have been carried out to demonstrate the performance of the proposed method; the corresponding simulation results show great improvement in image registration speed without losing repeatability. Finally, the CSA-SURF algorithm is applied to natural landmark registration task, showing that it works fine. The proposed method is a good candidate for image registration task for spacecraft autonomous navigation based on natural landmarks.



---

*Source: 8324298-2018-08-12.xml*