# Multiframe Astronomical Image Registration Based on Block Homography Estimation

Authors: Ran Li; Lin Luo; Yu Zhang
Journal: Journal of Sensors
 (2020)

---

## Abstract

Due to the influence of atmospheric turbulence, a time-variate video of an observed object by using the astronomical telescope drifts randomly with the passing of time. Thereafter, a series of images is obtained snapshotting from the video. In this paper, a method is proposed to improve the quality of astronomical images only through multiframe image registration and superimposition for the first time. In order to overcome the influence of anisoplanatism, a specific image registration algorithm based on multiple local homography transformations is proposed. Superimposing registered images can achieve an image with high definition. As a result, signal-to-noise ratio, contrast-to-noise ratio, and definition are improved significantly.


---

## Body

## 1. Introduction

Anisoplanatism continuously hinders the performance of spatial object observation by ground-based telescopes. Anisoplanatism is one of the key factors to directly influence spatial resolution in wide field imaging. In telescope imaging theory, the object spot keeps invariant in the scope of isoplanatic angle, and it has the same systematic point spread function (PSF), while if the imaging angle is out of the isoplanatic angle scope, the PSF is deemed to be varying with spatial position changes.In 1990, Labeyrie [1] published a paper to present a method on image restoration by spot inference. By using spot inference, Ozkan et al. [2] proposed a method to recover the extended spatial target phase spectra by using crossspectra; combining with the spot inference method, the extended astronomical image is restored. After that, Jin et al. [3] proposed a method to restore the extended spatial target phase by using bispectra, and the bispectra method is spatially invariant compared with the crossspectra method. Ayers and Dainty [4] innovatively proposed an image restoration concept by using blind deconvolution algorithm in the spatial target observation field. Zhang [5] adopted an iterative method proposed by Dainty to perform the calculations between spatial and spectral domains under certain restrictions and achieved a PSF that can map an image with the highest similarity to the real situation. Wang et al. [6] optimized Dainty’s method and added a filter to enhance robustness to noise, whereas in real spatial target observation procedure, a wide field extended spatial target is most likely to appear, and the visual field is much wider than the isoplanatic angle scope. From the experiment, the image restoration performance is impaired severely after the process of spot imaging or blind deconvolution. In spite of imaging via a self-adaptive optical telescope, the anisoplanatic effect is still salient [7, 8]. In order to overcome the shortages as mentioned above, Paxman et al. [9] proposed a method to simulate the anisoplanatic effect by using a wavefront coding technique. Thereafter, multiple frame superimposition was adopted instead of long-term exposure imaging to reduce the effect of anisoplanatism, and image definition was also improved [10]. Multiple local nonrigid transformations to impair spatial displacement and deformation among images were adopted [11, 12], as well as multiple images after transformation superimposition to realize image restoration. In order to improve image registration accuracy, establishing the image transformation model is necessary to proceed with various spatial transforms. Dan [13] proposed estimating the transformation model among short-exposure images with anisoplanatic effect based on neural networks, and the anisoplanatism is reduced by using image registration.The homography matrix is a kind of image spatial transformation model [14, 15]. Adopting homography transformation on each pixel of an image can realize linear or nonlinear transformation such as shifting, rotation, stretching, affine transformation, and prospect transformation [16, 17]. Some scholars propose using the block concept to calculate the homography matrix to achieve higher accuracy [18–20]. A coarse-to-fine method involves adopting image registration on the whole field to make the key points well-distributed. After that, a more accurate image registration on local field is used as a block process based on multiscale analysis to achieve better alignment performance [21]. In the procedure of local image registration, the block process proceeds based on the clustering method. A spectral clustering algorithm based on stratification consensus is proposed as an optimized clustering strategy to partition the image into several blocks [22]. A swarm optimization algorithm is proposed to adjust the clustering center position iteratively. Thus, the identification and mutual independence of each cluster are enhanced [23, 24]. The feature distribution is always overlooked in the procedure of block image registration. Deng et al. [25] propose a method that involves obtaining the feature distribution as a first step. Then, the image is partitioned into a number of blocks using clustering. Next, the feature points are assembled to the cluster centers along the margin directions and form a certain number of clusters. Lastly, block image registration on cluster scope between each two corresponding clusters is adopted to obtain a better image registration performance.This paper considers the multiframe image registration and superimposition concept to improve the image quality. Superimposing multiple low-resolution and noisy images can reduce the background noise and random drift. Image registration is deemed as an effective method to estimate displacements and deformations between two images. Thus, the transformation model is established after such estimation.Our contributions are as follows: (1) the effects of atmospheric turbulence cause astronomical images taken from the ground to have very little texture detail, and edge structures are not obvious and contain a lot of background noise. This makes it difficult for traditional algorithms (Scale Invariant Feature Transform (SIFT) [26], Speeded Up Robust Features (SURF) [27], Binary Robust Invariant Scalable Keypoints (BRISK) [28], and Features from Accelerated Segment Test (FAST) [29]) to extract feature points from astronomical images. In this paper, the Adaptive and Generic Accelerated Segment Test (AGAST) [30] algorithm and the DAISY [31] descriptor are used to extract and describe feature points of the image to overcome the above-mentioned shortcomings, so that a large number of key points with uniform distribution can be detected while maintaining high computing efficiency. (2) In order to further improve the calculation efficiency, this paper proposes a method based on the maximum entropy model for the selection of feature points and uses a method based on third-order statistics that is better than the Principal Components Analysis (PCA) algorithm [32, 33] to reduce the dimension of feature description vector. (3) The Fuzzy C-means (FCM) clustering is applied to feature points, and multiple block homography matrices are calculated between template image and image to be registered by the improved Random Sample Consensus (RANSAC) algorithm [34–36]. (4) The image is registered in blocks and superimposed to improve the quality of astronomical images.The remaining part of this paper is organized as follows. Section2 presents theories and discussion of the adopted method. Section 3 describes the experiment analysis, which analyzes the evaluation procedure of results processed by the proposed methods. Section 4 concludes with the strategies used in and theoretical basis of this paper.


## 2. Methods

The key issue in determining the relationship between two images captured at different times is establishing the transformation model. The geometrical transformation model is obtained via corresponding fixed points and their feature vectors or descriptors. Taking the unrobustness of feature point descriptors into consideration, an image matching error may occur. In this paper, the geometrical transformation model is calculated by homography estimation to solve the problem of matching errors. A point set is presented as(1)U=p1,p1′,p2,p2′,⋯,pi,pi′,where pi=xi,yi,1T and pi′=xi′,yi′,1T are two points in one image.Homography matrixH can be used to describe the spatial corresponding relationship between two images, as shown in Figure 1. The two points in two images can be transformed by each other via planeπ, presented as pi′=Hpi. The calculation of the homography matrix is done by selecting four points in each image, respectively, and a linear equation set is obtained. Then, the parameters of the homography matrix are obtained by soling such an equation set.Figure 1
The spatial corresponding relationship between two images.A proper homography transformation model is judged by measuring the Euclidean distance among the fixed points. The distance is deemed as a matching error, presented as(2)EHi2=dxi′,Hxi2+dxi,H−1xi′2,where d⋅ denotes the Euclidean distance between the fixed-point set and its corresponding set.In this paper, the traditional image registration algorithm based on the homography matrix is further improved. A multiframe astronomical image registration algorithm based on block homography estimation is proposed. Figure2 is the algorithm flow chart of the whole astronomical image registration process.Figure 2
Algorithm flow chart of the whole astronomical image registration process.


### 2.1. Reduced Feature Descriptor Calculation

The extraction and description procedures are preliminaries of homography estimation. This paper adopts an AGAST-DAISY algorithm on the template image and input image, respectively. With astronomical images, AGAST operator may detect a large number of key points. Therefore, a maximum entropy model-based method is used to refine the key points, and much of them will be reduced. Firstly, we take all the pixels in a whole image into consideration. In particular, a boundary piecewise process is used to deal with the pixels on the image boundaries, as shown in Figure3. Then, the conditional probability distribution model on the neighbor of each pixel is calculated, as Equation (3). Assume that each key point yields to Gaussian distribution on its neighbor.
(3)Px,y=12πσexp−fx,y−μ22σ2,where x,y denotes the position of one pixel. μ and σ2 denote the mean value and variation of one pixel on its neighbor. Assume that the mean value is approximate to the expected value. PX and Epf are the marginal probability and the expected value of axis X, respectively, shown as Equation (4) and Equation (5). The variation of Equation (3) can be presented as Equation (6).
(4)∑YPXPx,y=PX.Figure 3
The neighbor and the extended pixels of a boundary pixel under the piecewise concept.Solving outPX in Equation (4), the marginal probability distribution can be obtained. Putting the solved PX into Equation (5), the expected value is calculated as follows:
(5)Epf=∑x,yPx,yfx,y≈∑x,yPXPx,yfx,y,(6)σ2=∑x,yfx,y−Epf2.After the calculation of the conditional probability distribution of each key point, the entropy of one pixel on its neighbor is obtained as well, as shown in the following equation:(7)HP=−∑x,yPx,ylog10Px,y=−∑x,yPXPx,ylog10Px,y.A certain part of the key points with the highest entropy value is selected, and these pixels are deemed the refined key points, as shown in Figure4.Figure 4
(a) Satellite image; (b) AGAST key points of satellite image; (c) refined AGAST key points of satellite image; (d) rotated satellite image; (e) AGAST key points of rotated satellite image; (f) refined AGAST key points of rotated satellite image.
(a)(b)(c)(d)(e)(f)Next, we adopt the DAISY operator to generator feature descriptors for the key points. Each key point has its descriptor with the dimension of1×200. Thus, dimension reduction is necessary to be conducted for the consideration of efficiency. In order to overcome the shortcomings of the PCA algorithm on dimension reduction applications, we propose a method based on third-order statistics to conduct dimension reduction. To a feature vector with the dimension of 1×N, the three-order accumulation is calculated as shown in the following equation:
(8)c3xNm1,m2=1N∑n=1Nxnxn+m1xn+m2,where when n<1 or n>N, xn=0, and n is an integer number, m1,m2∈1,N. Then, the bispectrum of the feature vector is calculated by adopting Fourier transformation on the three-order accumulation. In particular, the frequency range of Fourier transform is set to 0,512 Hz, and Fourier transform is adopted on such frequency intervals, 0,64 Hz, 64,128 Hz,…, and 448,512 Hz. Thus, the feature vector bispectrum on 8 subfrequency intervals is calculated as
(9)Bxω1,ω2=∑m1=−∞+∞c3xNm1,m2e−jω1m1+ω2m2,where ω1 and ω2 are the spectral range of the bispectrum. The calculated bispectrum on 8 subfrequency bands is presented as B1ω1,ω2,B2ω1,ω2,⋯,B8ω1,ω2. Next, take the diagonal slice of each subfrequency band of the bispectrum, and select the maximum value of such diagonal slice as the element to form the feature vector. Thus, each feature vector contains 8 elements, shown as Equation (10) and Equation (11). Consequently, the dimension of the descriptor is reduced to 1×8,
(10)Siω=Biω1,ω2ω1=ω2,(11)Xi=argmaxSiω,where i=1,⋯,8 and Xi is the maximum value of the i′th subfrequency band Siω of the bispectrum. Therefore, the descriptor can be presented as X=X1,⋯,X8.From Figures5–7, the results obtained by using PCA have poor identity and discrimination among the existing features. Thus, it is difficult to distinguish certain features after the process by PCA, while the three-order accumulation processed feature vectors have much better identity and discrimination which can be used to figure out the certain features easily. Therefore, using three-order accumulation processed feature vectors in feature matching can achieve a superior performance.Figure 5
DAISY descriptor of one key point in satellite image: (a) the DAISY descriptors of two key points at the same position in two different satellite images; (b) the DAISY descriptors of two key points at the different positions in one satellite image.
(a)(b)Figure 6
Descriptor after PCA dimension reduction of one key point in satellite image: (a) the PCA processed descriptors of two key points at the same position in two different satellite images; (b) the PCA processed descriptors of two key points at the different positions in one satellite image.
(a)(b)Figure 7
DAISY descriptor proposed by our proposed three-order accumulation of one key point in satellite image: (a) the three-order accumulation processed descriptors of two key points at the same position in two different satellite images; (b) the three-order accumulation processed descriptors of two key points at the different positions in one satellite image.
(a)(b)


### 2.2. Image Registration Based on Block Homography Estimation

The determination of the number of clusters plays an important role in clustering. We use the Jeffrey divergence (JD) method to determine the number of clusters. The JD method is set as an object function in optimized polynomial form, shown as Equation (12). After recursive calculation of clustering, the JD value reaches its minimum value while the gap value reaches its maximum value simultaneously. The maximum intracluster compactness and the maximum intercluster discrimination are traded off, as well as ∂JD/∂n=0. The relationship curve between the gap value and number of clusters is shown in Figure 8(a). We record 100 iterations of gap calculation, and the maximum gap with respect to the number of clusters is shown in Figure 8(b).Figure 8
The relationship curve of gap statistics with respect to number of clusters of satellite image: (a) the relationship curve between the gap and the number of clusters; (b) the maximum gap with respect to the number of clusters; (c) the relationship curve of log of square summation of the number of points in one cluster with respect to the number of clusters.
(a)(b)(c)The relationship curve between the gap value and number of clusters is shown in Figure8(a). After 100 iterative calculations of gap, the result shows that the gap value reaches the maximum value when the number of clusters is 4, c=4, and the gap value is greater than 60. Thus, the most proper clustering can be performed under such conditions. As a result, the maximum intracluster compactness and the minimum intercluster overlap are achieved simultaneously, as shown in Figure 8(b). The log of square summation of the number of points in one cluster reaches the minimum value when c=4, which means the highest correlation among the intracluster points is achieved under such conditions, as shown in Figure 8(c):
(12)JDCi,Cj=12trΣi−1Σj+trΣj−1Σi+12μi−μjTΣi−1+Σjμi−μj,(13)Σi=∑u=1Nixu−μiTxu−μi,(14)μi=1Ni∑u=1Nixi,where i and j stand for the different clusters after clustering. Σi is the covariance matrix of the i′th cluster. μi is the mean value of the pixels in the i′th cluster.Determination of the number of clusters is a prerequisite of conducting clustering. In this paper, the number of target points is assumed to beN, presented as X1,X2,⋯,XN. The cluster centers are formed by using clustering, presented as C1,C2,⋯,Cj,⋯,Cc. In particular, the object function of clustering can be deemed to be the multiple target optimization procedure. In this paper, we set the FCM clustering function with the minimum entropy model as the object function. The optimized object function Jm is obtained when the two components reach the minimum value simultaneously, as shown in
(15)Jm=∑i=1N∑j=1CuijmXi−Cj2−1N∑i=1Nlog101C⋅VK∑j=1j≠iCΚXi−Cj,(16)uijk+1=1∑k=1CXi−Cjk/Xi−Ckk2/m−1,(17)Cjk=∑i=1Nuijkm⋅Xi∑i=1Nuijkm,(18)∑j=1Cuij=1,where k denotes the number of iterations, 1≤k<∞. uijis the membership function, which represents the probability of thei′th point,i=1,⋯,N, which belongs to thej′th cluster, as shown in Equation (16). Equation (17) represents the relationship between Xi and the j′th cluster center Cj, where Xi is the descriptor of the i′th key point with the dimension of 1×8. Cj denotes the descriptor of the j′th cluster center. N is the number of key points. Equations ((15)–(17)) are the iterative optimization procedures. Firstly, set the initial value of the membership function as uij0, which is a random value in the interval of 0,1. Then, set the terminate condition ε, ε∈0,1. When object function Jm fulfills the condition maxuijk+1−uijk2≤ε, the iteration is terminated, and the optimized object function Jm is obtained. Thereafter, all the cluster centers C1,C2,⋯,Cc are calculated, and the fuzzy clustering procedure is completed. The maximum correlation principle is adopted, as shown in Equation (18), on the clusters generated via FCM clustering. Then, the points with the maximum correlation to the cluster centers are found. In Figure 9, these found points are classified into the equivalent number of classes to the clusters. Lastly, blockwise homography estimation ρX,Y is conducted to realize local image registration, as shown in Equation (19), where X and Y represent the DAISY descriptor in two images, respectively.
(19)ρX,Y=∑XY−∑X∑Y/N∑X2−∑X2/N∑Y2−∑Y2/N.Figure 9
Key point clustering of satellite image: (a) in the left to right order, template satellite image, input satellite image, and two images with their clustered key points; (b) the corresponding key point matches of satellite image.
(a)(b)


### 2.3. Fixed Points for Homography Estimating

The homography matrix is a3×3 matrix with 8 degrees of freedom (DoF). Therefore, solving the 8 unknowns requires at least 4 noncolinear point positions in the coordinates. In this paper, it is necessary to conduct the homography estimation process as stated above on each cluster or block. Here, the neighbor around one pixel with the size of 3×3 is taken, and the probability distribution of this pixel on this neighbor area is calculated. Then, the entropy is obtained. The points with Nmax maximum entropy were selected. The value Nmax is calculated via
(20)CNmax4≥num_iter,where num_iter=2num_fixPoints×10 and num_fixPoints are the number of fixed points. Thus, Nmax needs to fulfill the condition Nmax≥num_iter. Here, num_fixPoints=4 andnum_iter=160. Select Nmax in the interval of 16~20. Thus, the iteration time of RANSAC is reduced from 160 to 20, and the computation efficiency is improved significantly. Equally, this optimization procedure is used on each corresponding cluster pair between two images, respectively. Four fixed points among theNmaxpoints in each cluster are selected, and the homography matrix is obtained via Equation (21) and Equation (22). Next, as shown in Figure 10, the point set X20 and X2′0 is calculated under the Levenberg-Marquardt principle to obtain the total error of forward and backward propagation, as shown in Equation (23):
(21)X2′0=H0X20,(22)X3′′0=H0−1X2′−1,(23)Disk=∑i=1Nx′ki−Hxik−12+xik−H−1x′ik−12.Figure 10
The summed propagation error between the template satellite image and the input satellite image: (a) the points in the first cluster; (b) the points in the second cluster; (c) the points in the third cluster; (d) the points in the fourth cluster.
(a)(b)(c)(d)The calculation with Equations ((21)–(23)) is performed recursively Nmax times. Therefore, the error between the point set and its mapping set at each time is presented as Dis0,Dis1,⋯,DisNmax, and the minimum value is selected, as shown in Equation (24). The calculated homography matrix is the optimized one via the above-stated method,
(24)Dis=minDis0,Dis1,⋯,DisN−1.Through the RANSAC algorithm, we can easily find out the best four feature points in each cluster to calculate the best homography matrix of this cluster. Through the geometric relationship of these four points, we can find out the center point of each cluster and then judge each pixel point in the whole image with these center points; the pixel is close to any cluster of center point, and the homography matrix is used to transform the pixel.


## 2.1. Reduced Feature Descriptor Calculation

The extraction and description procedures are preliminaries of homography estimation. This paper adopts an AGAST-DAISY algorithm on the template image and input image, respectively. With astronomical images, AGAST operator may detect a large number of key points. Therefore, a maximum entropy model-based method is used to refine the key points, and much of them will be reduced. Firstly, we take all the pixels in a whole image into consideration. In particular, a boundary piecewise process is used to deal with the pixels on the image boundaries, as shown in Figure3. Then, the conditional probability distribution model on the neighbor of each pixel is calculated, as Equation (3). Assume that each key point yields to Gaussian distribution on its neighbor.
(3)Px,y=12πσexp−fx,y−μ22σ2,where x,y denotes the position of one pixel. μ and σ2 denote the mean value and variation of one pixel on its neighbor. Assume that the mean value is approximate to the expected value. PX and Epf are the marginal probability and the expected value of axis X, respectively, shown as Equation (4) and Equation (5). The variation of Equation (3) can be presented as Equation (6).
(4)∑YPXPx,y=PX.Figure 3
The neighbor and the extended pixels of a boundary pixel under the piecewise concept.Solving outPX in Equation (4), the marginal probability distribution can be obtained. Putting the solved PX into Equation (5), the expected value is calculated as follows:
(5)Epf=∑x,yPx,yfx,y≈∑x,yPXPx,yfx,y,(6)σ2=∑x,yfx,y−Epf2.After the calculation of the conditional probability distribution of each key point, the entropy of one pixel on its neighbor is obtained as well, as shown in the following equation:(7)HP=−∑x,yPx,ylog10Px,y=−∑x,yPXPx,ylog10Px,y.A certain part of the key points with the highest entropy value is selected, and these pixels are deemed the refined key points, as shown in Figure4.Figure 4
(a) Satellite image; (b) AGAST key points of satellite image; (c) refined AGAST key points of satellite image; (d) rotated satellite image; (e) AGAST key points of rotated satellite image; (f) refined AGAST key points of rotated satellite image.
(a)(b)(c)(d)(e)(f)Next, we adopt the DAISY operator to generator feature descriptors for the key points. Each key point has its descriptor with the dimension of1×200. Thus, dimension reduction is necessary to be conducted for the consideration of efficiency. In order to overcome the shortcomings of the PCA algorithm on dimension reduction applications, we propose a method based on third-order statistics to conduct dimension reduction. To a feature vector with the dimension of 1×N, the three-order accumulation is calculated as shown in the following equation:
(8)c3xNm1,m2=1N∑n=1Nxnxn+m1xn+m2,where when n<1 or n>N, xn=0, and n is an integer number, m1,m2∈1,N. Then, the bispectrum of the feature vector is calculated by adopting Fourier transformation on the three-order accumulation. In particular, the frequency range of Fourier transform is set to 0,512 Hz, and Fourier transform is adopted on such frequency intervals, 0,64 Hz, 64,128 Hz,…, and 448,512 Hz. Thus, the feature vector bispectrum on 8 subfrequency intervals is calculated as
(9)Bxω1,ω2=∑m1=−∞+∞c3xNm1,m2e−jω1m1+ω2m2,where ω1 and ω2 are the spectral range of the bispectrum. The calculated bispectrum on 8 subfrequency bands is presented as B1ω1,ω2,B2ω1,ω2,⋯,B8ω1,ω2. Next, take the diagonal slice of each subfrequency band of the bispectrum, and select the maximum value of such diagonal slice as the element to form the feature vector. Thus, each feature vector contains 8 elements, shown as Equation (10) and Equation (11). Consequently, the dimension of the descriptor is reduced to 1×8,
(10)Siω=Biω1,ω2ω1=ω2,(11)Xi=argmaxSiω,where i=1,⋯,8 and Xi is the maximum value of the i′th subfrequency band Siω of the bispectrum. Therefore, the descriptor can be presented as X=X1,⋯,X8.From Figures5–7, the results obtained by using PCA have poor identity and discrimination among the existing features. Thus, it is difficult to distinguish certain features after the process by PCA, while the three-order accumulation processed feature vectors have much better identity and discrimination which can be used to figure out the certain features easily. Therefore, using three-order accumulation processed feature vectors in feature matching can achieve a superior performance.Figure 5
DAISY descriptor of one key point in satellite image: (a) the DAISY descriptors of two key points at the same position in two different satellite images; (b) the DAISY descriptors of two key points at the different positions in one satellite image.
(a)(b)Figure 6
Descriptor after PCA dimension reduction of one key point in satellite image: (a) the PCA processed descriptors of two key points at the same position in two different satellite images; (b) the PCA processed descriptors of two key points at the different positions in one satellite image.
(a)(b)Figure 7
DAISY descriptor proposed by our proposed three-order accumulation of one key point in satellite image: (a) the three-order accumulation processed descriptors of two key points at the same position in two different satellite images; (b) the three-order accumulation processed descriptors of two key points at the different positions in one satellite image.
(a)(b)


## 2.2. Image Registration Based on Block Homography Estimation

The determination of the number of clusters plays an important role in clustering. We use the Jeffrey divergence (JD) method to determine the number of clusters. The JD method is set as an object function in optimized polynomial form, shown as Equation (12). After recursive calculation of clustering, the JD value reaches its minimum value while the gap value reaches its maximum value simultaneously. The maximum intracluster compactness and the maximum intercluster discrimination are traded off, as well as ∂JD/∂n=0. The relationship curve between the gap value and number of clusters is shown in Figure 8(a). We record 100 iterations of gap calculation, and the maximum gap with respect to the number of clusters is shown in Figure 8(b).Figure 8
The relationship curve of gap statistics with respect to number of clusters of satellite image: (a) the relationship curve between the gap and the number of clusters; (b) the maximum gap with respect to the number of clusters; (c) the relationship curve of log of square summation of the number of points in one cluster with respect to the number of clusters.
(a)(b)(c)The relationship curve between the gap value and number of clusters is shown in Figure8(a). After 100 iterative calculations of gap, the result shows that the gap value reaches the maximum value when the number of clusters is 4, c=4, and the gap value is greater than 60. Thus, the most proper clustering can be performed under such conditions. As a result, the maximum intracluster compactness and the minimum intercluster overlap are achieved simultaneously, as shown in Figure 8(b). The log of square summation of the number of points in one cluster reaches the minimum value when c=4, which means the highest correlation among the intracluster points is achieved under such conditions, as shown in Figure 8(c):
(12)JDCi,Cj=12trΣi−1Σj+trΣj−1Σi+12μi−μjTΣi−1+Σjμi−μj,(13)Σi=∑u=1Nixu−μiTxu−μi,(14)μi=1Ni∑u=1Nixi,where i and j stand for the different clusters after clustering. Σi is the covariance matrix of the i′th cluster. μi is the mean value of the pixels in the i′th cluster.Determination of the number of clusters is a prerequisite of conducting clustering. In this paper, the number of target points is assumed to beN, presented as X1,X2,⋯,XN. The cluster centers are formed by using clustering, presented as C1,C2,⋯,Cj,⋯,Cc. In particular, the object function of clustering can be deemed to be the multiple target optimization procedure. In this paper, we set the FCM clustering function with the minimum entropy model as the object function. The optimized object function Jm is obtained when the two components reach the minimum value simultaneously, as shown in
(15)Jm=∑i=1N∑j=1CuijmXi−Cj2−1N∑i=1Nlog101C⋅VK∑j=1j≠iCΚXi−Cj,(16)uijk+1=1∑k=1CXi−Cjk/Xi−Ckk2/m−1,(17)Cjk=∑i=1Nuijkm⋅Xi∑i=1Nuijkm,(18)∑j=1Cuij=1,where k denotes the number of iterations, 1≤k<∞. uijis the membership function, which represents the probability of thei′th point,i=1,⋯,N, which belongs to thej′th cluster, as shown in Equation (16). Equation (17) represents the relationship between Xi and the j′th cluster center Cj, where Xi is the descriptor of the i′th key point with the dimension of 1×8. Cj denotes the descriptor of the j′th cluster center. N is the number of key points. Equations ((15)–(17)) are the iterative optimization procedures. Firstly, set the initial value of the membership function as uij0, which is a random value in the interval of 0,1. Then, set the terminate condition ε, ε∈0,1. When object function Jm fulfills the condition maxuijk+1−uijk2≤ε, the iteration is terminated, and the optimized object function Jm is obtained. Thereafter, all the cluster centers C1,C2,⋯,Cc are calculated, and the fuzzy clustering procedure is completed. The maximum correlation principle is adopted, as shown in Equation (18), on the clusters generated via FCM clustering. Then, the points with the maximum correlation to the cluster centers are found. In Figure 9, these found points are classified into the equivalent number of classes to the clusters. Lastly, blockwise homography estimation ρX,Y is conducted to realize local image registration, as shown in Equation (19), where X and Y represent the DAISY descriptor in two images, respectively.
(19)ρX,Y=∑XY−∑X∑Y/N∑X2−∑X2/N∑Y2−∑Y2/N.Figure 9
Key point clustering of satellite image: (a) in the left to right order, template satellite image, input satellite image, and two images with their clustered key points; (b) the corresponding key point matches of satellite image.
(a)(b)


## 2.3. Fixed Points for Homography Estimating

The homography matrix is a3×3 matrix with 8 degrees of freedom (DoF). Therefore, solving the 8 unknowns requires at least 4 noncolinear point positions in the coordinates. In this paper, it is necessary to conduct the homography estimation process as stated above on each cluster or block. Here, the neighbor around one pixel with the size of 3×3 is taken, and the probability distribution of this pixel on this neighbor area is calculated. Then, the entropy is obtained. The points with Nmax maximum entropy were selected. The value Nmax is calculated via
(20)CNmax4≥num_iter,where num_iter=2num_fixPoints×10 and num_fixPoints are the number of fixed points. Thus, Nmax needs to fulfill the condition Nmax≥num_iter. Here, num_fixPoints=4 andnum_iter=160. Select Nmax in the interval of 16~20. Thus, the iteration time of RANSAC is reduced from 160 to 20, and the computation efficiency is improved significantly. Equally, this optimization procedure is used on each corresponding cluster pair between two images, respectively. Four fixed points among theNmaxpoints in each cluster are selected, and the homography matrix is obtained via Equation (21) and Equation (22). Next, as shown in Figure 10, the point set X20 and X2′0 is calculated under the Levenberg-Marquardt principle to obtain the total error of forward and backward propagation, as shown in Equation (23):
(21)X2′0=H0X20,(22)X3′′0=H0−1X2′−1,(23)Disk=∑i=1Nx′ki−Hxik−12+xik−H−1x′ik−12.Figure 10
The summed propagation error between the template satellite image and the input satellite image: (a) the points in the first cluster; (b) the points in the second cluster; (c) the points in the third cluster; (d) the points in the fourth cluster.
(a)(b)(c)(d)The calculation with Equations ((21)–(23)) is performed recursively Nmax times. Therefore, the error between the point set and its mapping set at each time is presented as Dis0,Dis1,⋯,DisNmax, and the minimum value is selected, as shown in Equation (24). The calculated homography matrix is the optimized one via the above-stated method,
(24)Dis=minDis0,Dis1,⋯,DisN−1.Through the RANSAC algorithm, we can easily find out the best four feature points in each cluster to calculate the best homography matrix of this cluster. Through the geometric relationship of these four points, we can find out the center point of each cluster and then judge each pixel point in the whole image with these center points; the pixel is close to any cluster of center point, and the homography matrix is used to transform the pixel.


## 3. Experimental Results and Analysis

### 3.1. Data Acquisition

In the experiment, we captured 10 frames of Saturn images and Mars images by snapshotting from the video with the same time interval. Due to the effect of atmospheric turbulence, the snapshotted images are blurred intensively. Meanwhile, sophisticated deformations exist in the images as well. The imaging system is set up as follows: the lens diameter of the telescope isD=1.8 m; atmospheric turbulence is presented as λl′=9.52×10−6, where l′=13.6 m. The wavelength of the visible light is λ=0.7×10−7. Thus, the cutoff frequency isfc=D/λl′=1.89×105 Hz. The size of the captured Saturn images is 1024×1024, while the size of the Mars images is 512×512.


### 3.2. Evaluation Standard

In this paper, signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), root mean square error (RMSE), and power spectrum are adopted to illustrate the performance of using different methods. As shown in Equations ((25)–(29)),
(25)SNR=10×log10S2N2,where S denotes the useful signal. N denotes the noise in one image,
(26)CNR=S¯t−S¯bσb,where S¯t represents the mean pixel value of the target area. S¯b represents the mean pixel value of the background area. σb is the standard deviation of the pixels in background area,
(27)RMSE=∑im∑jnfi,j−f¯2m×n,where m and n represent the size of an image. f¯ is the mean value of the pixels of a whole image.The autocorrelation matrix of an imagefx,y is shown as
(28)Rffτ,ς=∑x=1M∑y=1Nfτ,ς⋅fx−τ,y−ς.The logarithm power spectrum is presented as(29)log10Ef2u,v=Rff0,0=12π∫−∞+∞∫−∞+∞Sffjω,jυdωdυ,where fu,v stands for the pixel value of the point at the position u,v in the coordinates. Sffjω,jυ is the matrix after the Fourier transformation on the autocorrelation matrixRffτ,ς.In this paper, we adopt the AGAST operator to detect key points in Saturn images and Mars images. Then, the key points are refined by using the maximum entropy method, and the refining ratio is set to 0.3 to make all the points located at the most feature-representative area, as shown in Figure11. In the order of left to right of each row are the raw data, the image with AGAST key points, and the image with refined AGAST key points. This paper adopts a block homography estimation method based on FCM clustering. The images processed by FCM clustering form a certain number of clusters, as shown in Figure 12. In Figures 12(a) and 12(b), from left to right are the raw data and the key points in the first image and second image after clustering. Different colors denote different clusters. Each class of the key points forms a certain block.Figure 11
(a) Saturn image; (b) AGAST key points in Saturn image; (c) refined AGAST key points in Saturn image; (d) Mars image; (e) AGAST key points in Mars image; (f) refined AGAST key points in Mars image.
(a)(b)(c)(d)(e)(f)Figure 12
Astronomical image key point clustering. (a) In the order of left to right are the raw data, clustered key points in one Saturn image, and clustered key points in another Saturn image. (b) In the order of left to right are the raw data, clustered key points in one Mars image, and clustered key points in another Mars image.
(a)(b)


### 3.3. Evaluation Standard of Noise

Figures13 and 14 are the processed results of Saturn image and Mars image, respectively, and show two images captured at different times, 10 image superimposition, 3 images processed by the proposed method superimposition, 5 images processed by the proposed method superimposition, and 10 images processed by the proposed method superimposition. Figures 13(g) and 14(g) show the zoomed in images of the part framed by the green box in Figures 13(a)–13(f) and Figures 14(a)–14(f). It can be seen that the image noise is being depressed with the increasing of the superimposed number of images.Figure 13
The image registration result of Saturn images: (a) original image; (b) another original image; (c) 10 original image superimposition (d) 3 images processed by the proposed method superimposition; (e) 5 images processed by proposed method superimposition; (f) 10 images processed by the proposed method superimposition; (g) the zoomed in images of the part framed by green box in (a–f).
(a)(b)(c)(d)(e)(f)(g)Figure 14
The image registration result of Mars images: (a) original image; (b) another original image; (c) 10 original image superimposition; (d) 3 images processed by the proposed method superimposition; (e) 5 images processed by the proposed method superimposition; (f) 10 images processed by the proposed method superimposition; (g) the zoomed in images of the part framed by green box in (a–f).
(a)(b)(c)(d)(e)(f)(g)Figures15 and 16 show that, with the increasing of the number of the superimposed images, SNR and CNR are increasing, while the RMSE is decreasing. However, the evaluation index tends to become steady when the superimposed number of images reaches 5, as shown in Tables 1 and 2.Figure 15
The relationship curve of evaluation index with respect to the number of superimposed images: (a) the relationship curve of SNR with respect to the number of superimposed images; (b) the relationship curve of CNR with respect to the number of superimposed images; (c) the relationship curve of RMSE with respect to the number of superimposed images.
(a)(b)(c)Figure 16
The relationship curve of evaluation index with respect to the number of superimposed images: (a) the relationship curve of SNR with respect to the number of superimposed images; (b) the relationship curve of CNR with respect to the number of superimposed images; (c) the relationship curve of RMSE with respect to the number of superimposed images.
(a)(b)(c)Table 1
All the evaluation indexes for Saturn image.
Evaluation indexesSNR (dB)CNRRMSERunning time (s)Saturn images processed by the proposed method superimposition3 images21.067614.27630.025515.38615 images22.072116.03820.018024.934710 images22.298816.46320.017246.5839Original Saturn images superimposed3 images14.63357.13230.0425—5 images16.72607.59950.0369—10 images16.85617.54530.0367—Table 2
All the evaluation indexes for Mars image.
Evaluation indexesSNR (dB)CNRRMSERunning time (s)Mars images processed by the proposed method superimposition3 images20.0721103.00300.02833.76415 images20.4280112.07030.02706.483710 images20.5535113.86790.026311.8962Original Mars images superimposed3 images14.772370.55780.0510—5 images15.167274.36500.0473—10 images15.515571.72830.0494—


### 3.4. Power Spectrum Analysis

Figures17–20 show the logarithm power spectrum of the Saturn image, one-dimension logarithm power spectrum of the Saturn image, the logarithm power spectrum of the Mars image, and the one-dimension logarithm power spectrum of the Mars image. As can be seen in Figures 18 and 20, the intensity of the one-dimension logarithm power spectrum of the image processed by our proposed methods is lower than that of the original image and the image processed by some classic methods. Therefore, more energy is centered on the low-frequency domain. Consequently, the proposed algorithm outperforms the other methods even the state of the art on image quality enhancement.Figure 17
The power spectrum of the Saturn image: (a) original image; (b) 10 original image superimposition; (c) 3 images processed by the proposed method superimposition; (d) 5 images processed by the proposed method superimposition; (e) 10 images processed by the proposed method superimposition.
(a)(b)(c)(d)(e)Figure 18
The one-dimension logarithm power spectrum of the Saturn image. (a) The red, green, and blue curves represent the one-dimension logarithm power spectrum of original image, 10 original image superimposition, and 10 images processed by the proposed method superimposition, respectively; (b) the red, green, and blue curves represent the one-dimension logarithm power spectrum of 3 images, 5 images, and 10 images processed by the proposed method superimposition.
(a)(b)Figure 19
The power spectrum of the Mars image: (a) original image; (b) 10 original image superimposition; (c) 3 images processed by the proposed method superimposition; (d) 5 images processed by the proposed method superimposition; (e) 10 images processed by the proposed method superimposition.
(a)(b)(c)(d)(e)Figure 20
The one-dimension logarithm power spectrum of the Mars image. (a) The red, green, and blue curves represent the one-dimension logarithm power spectrum of original image, 10 original image superimposition, and 10 images processed by the proposed method superimposition, respectively; (b) the red, green, and blue curves represent the one-dimension logarithm power spectrum of 3 images, 5 images, and 10 images processed by the proposed method superimposition.
(a)(b)


## 3.1. Data Acquisition

In the experiment, we captured 10 frames of Saturn images and Mars images by snapshotting from the video with the same time interval. Due to the effect of atmospheric turbulence, the snapshotted images are blurred intensively. Meanwhile, sophisticated deformations exist in the images as well. The imaging system is set up as follows: the lens diameter of the telescope isD=1.8 m; atmospheric turbulence is presented as λl′=9.52×10−6, where l′=13.6 m. The wavelength of the visible light is λ=0.7×10−7. Thus, the cutoff frequency isfc=D/λl′=1.89×105 Hz. The size of the captured Saturn images is 1024×1024, while the size of the Mars images is 512×512.


## 3.2. Evaluation Standard

In this paper, signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), root mean square error (RMSE), and power spectrum are adopted to illustrate the performance of using different methods. As shown in Equations ((25)–(29)),
(25)SNR=10×log10S2N2,where S denotes the useful signal. N denotes the noise in one image,
(26)CNR=S¯t−S¯bσb,where S¯t represents the mean pixel value of the target area. S¯b represents the mean pixel value of the background area. σb is the standard deviation of the pixels in background area,
(27)RMSE=∑im∑jnfi,j−f¯2m×n,where m and n represent the size of an image. f¯ is the mean value of the pixels of a whole image.The autocorrelation matrix of an imagefx,y is shown as
(28)Rffτ,ς=∑x=1M∑y=1Nfτ,ς⋅fx−τ,y−ς.The logarithm power spectrum is presented as(29)log10Ef2u,v=Rff0,0=12π∫−∞+∞∫−∞+∞Sffjω,jυdωdυ,where fu,v stands for the pixel value of the point at the position u,v in the coordinates. Sffjω,jυ is the matrix after the Fourier transformation on the autocorrelation matrixRffτ,ς.In this paper, we adopt the AGAST operator to detect key points in Saturn images and Mars images. Then, the key points are refined by using the maximum entropy method, and the refining ratio is set to 0.3 to make all the points located at the most feature-representative area, as shown in Figure11. In the order of left to right of each row are the raw data, the image with AGAST key points, and the image with refined AGAST key points. This paper adopts a block homography estimation method based on FCM clustering. The images processed by FCM clustering form a certain number of clusters, as shown in Figure 12. In Figures 12(a) and 12(b), from left to right are the raw data and the key points in the first image and second image after clustering. Different colors denote different clusters. Each class of the key points forms a certain block.Figure 11
(a) Saturn image; (b) AGAST key points in Saturn image; (c) refined AGAST key points in Saturn image; (d) Mars image; (e) AGAST key points in Mars image; (f) refined AGAST key points in Mars image.
(a)(b)(c)(d)(e)(f)Figure 12
Astronomical image key point clustering. (a) In the order of left to right are the raw data, clustered key points in one Saturn image, and clustered key points in another Saturn image. (b) In the order of left to right are the raw data, clustered key points in one Mars image, and clustered key points in another Mars image.
(a)(b)


## 3.3. Evaluation Standard of Noise

Figures13 and 14 are the processed results of Saturn image and Mars image, respectively, and show two images captured at different times, 10 image superimposition, 3 images processed by the proposed method superimposition, 5 images processed by the proposed method superimposition, and 10 images processed by the proposed method superimposition. Figures 13(g) and 14(g) show the zoomed in images of the part framed by the green box in Figures 13(a)–13(f) and Figures 14(a)–14(f). It can be seen that the image noise is being depressed with the increasing of the superimposed number of images.Figure 13
The image registration result of Saturn images: (a) original image; (b) another original image; (c) 10 original image superimposition (d) 3 images processed by the proposed method superimposition; (e) 5 images processed by proposed method superimposition; (f) 10 images processed by the proposed method superimposition; (g) the zoomed in images of the part framed by green box in (a–f).
(a)(b)(c)(d)(e)(f)(g)Figure 14
The image registration result of Mars images: (a) original image; (b) another original image; (c) 10 original image superimposition; (d) 3 images processed by the proposed method superimposition; (e) 5 images processed by the proposed method superimposition; (f) 10 images processed by the proposed method superimposition; (g) the zoomed in images of the part framed by green box in (a–f).
(a)(b)(c)(d)(e)(f)(g)Figures15 and 16 show that, with the increasing of the number of the superimposed images, SNR and CNR are increasing, while the RMSE is decreasing. However, the evaluation index tends to become steady when the superimposed number of images reaches 5, as shown in Tables 1 and 2.Figure 15
The relationship curve of evaluation index with respect to the number of superimposed images: (a) the relationship curve of SNR with respect to the number of superimposed images; (b) the relationship curve of CNR with respect to the number of superimposed images; (c) the relationship curve of RMSE with respect to the number of superimposed images.
(a)(b)(c)Figure 16
The relationship curve of evaluation index with respect to the number of superimposed images: (a) the relationship curve of SNR with respect to the number of superimposed images; (b) the relationship curve of CNR with respect to the number of superimposed images; (c) the relationship curve of RMSE with respect to the number of superimposed images.
(a)(b)(c)Table 1
All the evaluation indexes for Saturn image.
Evaluation indexesSNR (dB)CNRRMSERunning time (s)Saturn images processed by the proposed method superimposition3 images21.067614.27630.025515.38615 images22.072116.03820.018024.934710 images22.298816.46320.017246.5839Original Saturn images superimposed3 images14.63357.13230.0425—5 images16.72607.59950.0369—10 images16.85617.54530.0367—Table 2
All the evaluation indexes for Mars image.
Evaluation indexesSNR (dB)CNRRMSERunning time (s)Mars images processed by the proposed method superimposition3 images20.0721103.00300.02833.76415 images20.4280112.07030.02706.483710 images20.5535113.86790.026311.8962Original Mars images superimposed3 images14.772370.55780.0510—5 images15.167274.36500.0473—10 images15.515571.72830.0494—


## 3.4. Power Spectrum Analysis

Figures17–20 show the logarithm power spectrum of the Saturn image, one-dimension logarithm power spectrum of the Saturn image, the logarithm power spectrum of the Mars image, and the one-dimension logarithm power spectrum of the Mars image. As can be seen in Figures 18 and 20, the intensity of the one-dimension logarithm power spectrum of the image processed by our proposed methods is lower than that of the original image and the image processed by some classic methods. Therefore, more energy is centered on the low-frequency domain. Consequently, the proposed algorithm outperforms the other methods even the state of the art on image quality enhancement.Figure 17
The power spectrum of the Saturn image: (a) original image; (b) 10 original image superimposition; (c) 3 images processed by the proposed method superimposition; (d) 5 images processed by the proposed method superimposition; (e) 10 images processed by the proposed method superimposition.
(a)(b)(c)(d)(e)Figure 18
The one-dimension logarithm power spectrum of the Saturn image. (a) The red, green, and blue curves represent the one-dimension logarithm power spectrum of original image, 10 original image superimposition, and 10 images processed by the proposed method superimposition, respectively; (b) the red, green, and blue curves represent the one-dimension logarithm power spectrum of 3 images, 5 images, and 10 images processed by the proposed method superimposition.
(a)(b)Figure 19
The power spectrum of the Mars image: (a) original image; (b) 10 original image superimposition; (c) 3 images processed by the proposed method superimposition; (d) 5 images processed by the proposed method superimposition; (e) 10 images processed by the proposed method superimposition.
(a)(b)(c)(d)(e)Figure 20
The one-dimension logarithm power spectrum of the Mars image. (a) The red, green, and blue curves represent the one-dimension logarithm power spectrum of original image, 10 original image superimposition, and 10 images processed by the proposed method superimposition, respectively; (b) the red, green, and blue curves represent the one-dimension logarithm power spectrum of 3 images, 5 images, and 10 images processed by the proposed method superimposition.
(a)(b)


## 4. Conclusion

In this paper, an image quality enhancement method based on multiple image registration and superimposition is analyzed and discussed. A series of images captured from a video captured via an astronomical telescope are aligned at different times. The major contribution of this paper is improving image registration accuracy and computation efficiency. As regards computation efficiency, the AGAST operator is first adopted to extract key points, and a key point refining method based on maximum entropy is used to reduce the large number of key points to a much smaller one. Then, we design a method by using the principle of high-order statistics to reduce the DAISY descriptors of densely located key points. Last, but not the least, a maximum entropy model-based method with RANSAC is used to decrease the iteration times in the homography estimation.To improve accuracy, a block process concept is adopted in this paper to conduct the multiple local image registration. The blocks are generated by FCM clustering. Specifically, the object function is set under the condition of minimum propagation error. The blocks are obtained via iterative calculation of cluster centers and maximum correlation to the cluster centers. As a result, the highest compactness of the intrablock points and relatively high independence of the interblock points are achieved. Moreover, multiple registered image superimposition has an effective impact on noise suppression and pixel drifting calibration, which substantially improves SNR and CNR.



---

*Source: 8849552-2020-12-10.xml*