A Real-Time Star Tailing Removal Method Based on Fast Blur Kernel Estimations

-e number of star points and the accuracy of star centroid extraction are the key factors that affect the performance of the star sensor under high dynamic conditions.-emotion blur results in the star point trailing, which consequently leads to the decline of star centroid extraction accuracy and in some cases lead to extraction failure. In order to improve the dynamic performance of the star sensor, in this work, we propose a real-time star trailing removal method based on fast blur kernel estimation. First, in order to minimize the influence of noise on parameter estimation, we use principal component analysis (PCA) in the dual-frequency spectrum domain to estimate the angle of blur kernel. In addition, an adjustable weighting method is proposed to estimate the length of blur kernel. So, we are able to quickly estimate the high-precision blur kernel based on a single degraded image. Moreover, an area filtering method based on the hyper-Laplacian prior recovery algorithm is also proposed. -is algorithm quickly reconstructs the star points in the tracking windows and effectively removes the star point tailing in real time. -e computational efficiency of the proposed algorithm is 5 times superior to the traditional method and 15 times superior to the existing accelerated iterative method.-e experimental results show that the proposed algorithm removes the trailing star quickly and effectively, under low SNR. In addition, the proposed method effectively improves the number of extracted star points and the accuracies of star centroids.


Introduction
Star sensor device forms the basis of celestial navigation. e three-axis attitude information of spacecraft is calculated by using the observed star's position [1]. e star sensor has three working modes: whole sky mode, local sky mode, and tracking mode. e tracking mode is the main mode of star sensor's working. e star sensor performs poorly under high dynamic conditions. is is due to the fast motion of the detector relative to the star. When the camera is capturing images during the fast movement of the vehicle, the resulting images are blurred. As presented in Figure 1, the motion blur disperses the energy of stars in the image, consequently, leading to the reduction in the signal-to-noise ratio (SNR).
is greatly reduces the star's centroid extraction accuracy and even fails to extract star points due to the noise submergence, thus resulting in tracking failure. e aforementioned factors determine if the star sensor is able to perform normally under high dynamic conditions. erefore, in order to improve the performance of the star sensor under dynamic conditions, it is essential to improve the number of star points and the accuracy of centroid extraction.
ere are two major techniques that are used to enhance the accuracy of star centroid extraction, namely, hardware enhancement and software restoration. e major purpose of hardware enhancement is to reduce the star blur by improving the circuit design. In [2], Bezooijen et al. proposed a delay integration method. is technique has the ability to suppress the motion blur by improving the imaging circuit of the star sensor to a certain extent. In addition, this technique also improves the image SNR. However, the circuit design for the implementation of this technique is very complex and results in the reduction of the imaging area of the sensor. Moreover, this technique is only able to cope with the motion blur in the vertical direction. us, we need additional software processing for achieving our goal. e software restoration method mainly removes the tail of the star spot for improving the accuracy of star's centroid extraction. e software restoration method is divided into two steps: blur kernel estimation and image reconstruction. In the blind restoration technique for blurred images, most of the methods use Radon transform [3][4][5][6], Hough transform, and other line detection techniques to estimate the motion blur angle in the frequency domain. e length estimation usually uses curve fitting [3], bispectrum [7], power spectrum [8], and cepstrum [9]. Few other methods also make use of bilinear interpolation to refine the blur length estimation [10]. At present, most of the efficient methods [11,12] to enhance the accuracy of star centroid extraction are based on the improvements made in Radon transform. e technique of Radon transform traverses all the angles on the spectrum chart to find the maximum value.
is maximum value is the estimated angle, as presented in Figure 2. e Radon transform is not only affected easily by image resolution and noise, but it also has poor robustness and higher computation cost.
In recent years, deep learning [13][14][15] has also begun to be applied to the parameter estimation of the blur kernel. In reference to the uncertainty of the actual degraded image blur type, the literature [16] uses convolutional neural networks (CNN) to map the input samples to the discriminative feature space in a supervised learning manner to realize the analysis and classification of the cause of the blur and call the definition regression neural network to efficiently estimate the blur kernel. Literature [17] uses CNN to mine motion blur features and uses support vector machine (SVM) to map them into blur angles and blur lengths to estimate blur kernel. Some algorithms [7,12,18] aim at nonvertical and nonhorizontal blurred images and simplify the estimation process through data dimensionality reduction, but it is easy to cause interpolation and data loss [19]. ere are also many research studies on restoration methods based on deep learning. e literature [20] combines the advantages of energy functional optimization and deep learning inverse convolution models and uses halfduadratic splitting (HQS) to decompose the energy functional and applies CNN noise prior deals with the subproblem of denoising. e Chinese University of Hong Kong researched the "from coarse to fine" image restoration process [21] and proposed a scale recurrent network that introduces coding-decoding residual network blocks at each scale. is network has fewer parameters and the training process is simplified accordingly. Literature [22] uses generative adversarial network (CAN) based on conditional CAN and multicomponent loss function. And using the game idea of generating confrontation makes the generator and the discriminator reach the Nash equilibrium, so as to restore the image. Although deep learning develops rapidly and has excellent results, it relies on a huge dataset and requires high hardware, so it is not suitable for star sensors with small computing power.
Most of the existing star tracking algorithms perform nonblind restoration.
is technique relies on auxiliary tools, such as the gyroscope of the star sensor to estimate the blur kernel to restore the star map for improving the accuracy of star centroid [23][24][25][26]. In [27], the authors propose a star centroid motion model according to the angular velocity obtained from the gyroscope. is estimation is then used to obtain the blur kernel and use Richardson Lucy (RL) algorithm [28,29] to restore star map. However, the accuracy of the gyroscope gradually declines with time. us, the long-term orbit operations by star sensors cannot be implemented. On the other hand, few studies suggest the direct estimation of blur kernel according to the star map to restore the star map. In [30], Jiang et al. proposed an improved Radon transform for the estimation of blur kernel. e method applies Z-function and double threshold mask before estimating the parameters to improve the accuracy of blur kernel estimation. However, the preprocessing of the proposed algorithm is computationally inefficient. is is because this method requires the complete image of the star to be processed, while the star sensor is in the tracking mode.
us, the algorithm is unable to meet the real-time tracking requirements of the star sensor. e process of reconstructing the image is to use inverse convolution to reverse the image degradation process. e energy functional inverse convolution model depends on the image prior. e difference between the different methods is roughly that the norm used in the model is different. e commonly used Wiener filter is through the L2 norm to constrain the image gradient to obey the Gaussian distribution, and the L2 norm is a convex function, so the Wiener filter is solved quickly, but it is prone to ringing. Literature [31] uses the L1 norm to constrain the image gradient to obey the Laplace distribution. e well-known ROF total variation (TV) denoising model [32] is also a typical application of this type of prior. Although the total variational norm is a convex function, it is not derivable. erefore, in recent years, its total variational inverse convolution model mainly focuses on how to optimize the solution. When the norm pε[0, 1], the constrained image gradient obeys the hyper-Laplace distribution, which is the prior closest to the heavy-tailed [33] distribution of natural images, especially when pϵ[0.5, 0.8] performs best. e commonly used methods of star image reconstruction include the Wiener filtering, RL algorithm, and regularization. In [26], Wang proposed an application of Wiener filter; however, the ringing effect in the proposed method is significant. Similarly, in [24], Wu presents the application of the regularization method; however, this method is also affected by the ringing effect. In [27], Sun makes an effort to suppress the ringing effect by using the RL reconstruction algorithm. However, the algorithm is computationally inefficient and takes a lot of time to execute. As time proceeds, the ringing effect increases as the number of iterations increases. Jiang in [30] uses second-order extrapolation to accelerate the RL method; however, the iteration process of the RL method still needs to be set manually.
According to the star image estimation method, the difficulty lies in the low accuracy of the blur kernel estimation and the poor real-time performance of the whole algorithm. As a result, it is difficult to remove the trailing star points quickly and effectively. erefore, it is necessary to implement star image estimation algorithm with high accuracy and good real-time performance. In this work, the principal component analysis (PCA) [34][35][36] is used in the dual-frequency spectral domain and an adjustable weighting method is proposed to estimate the blur kernel parameters of degraded images. A high-precision blur kernel is quickly estimated by using a single frame blur star map. Similarly, an area filtering method based on hyper-Laplacian prior [37] recovery algorithm is proposed. In this method, we are able to quickly reconstruct the star points in the tracking windows. e computation speed of the proposed algorithm is 5 times higher than that of the traditional method and 15 times higher than that of the existing accelerated iterative method proposed in [30]. e rest of the manuscript is organized as follows: in Section 2, we present the formation mechanism of star trailing. In Section 3, we present the proposed method for reconstructing the star points in tracking windows. In Section 4, we present the simulation results and analysis of the proposed method. Finally, in Section 5, we conclude our work.

Mechanism Analysis
e relative motion of the target star and the imaging device mounted on the star sensor in the limited exposure time is the source of the motion blur of the captured star images. In this section, we present the mathematical model for the motion blur embedded in the acquired star images. At rest, the energy distribution of stars in the image generally follows the Gaussian distribution E (x, y).
where σ represents the variance of energy distribution (σ is dependent on the parameters of star sensor) and E 0 represents the total energy of the stars. (x 0 , y 0 ) presents the center point. e star image f (x, y) is represented by When the star sensor operates in a high dynamic environment, there is a relative motion between the star and e star is a point light source at infinity. Considering the short exposure time and small rotation angle in the exposure time, the movement of the star spots in the image plane is approximated by linear motion [9]. e energy distribution E (x, y) of the star's image disperses on the sensor, which is in a relative motion to the star. is phenomenon leads to the formation of "tail" of the star spot. In addition, the center of the star spot is also displaced. e length L of motion blur is directly proportional to exposure time δt and speed v of the star sensor. Now, the energy distribution of star is where ∆x � x − x 0 − rcosα and ∆y � y − y 0 -rsinα. α, r, and L denote the motion angle, displacement, and star point motion length in the exposure time δt, respectively. e blurred image is represented as g(x, y) (g(x, y) � C·E b ). After the addition of additive noise n (x, y), the degraded image g(x, y) is expressed in the form of convolution operation as follows: where h (x, y) represents the motion blur kernel. is kernel contains the motion information and interacts with the noise in the image to get the degraded image. Figure 3 shows the block diagram for the proposed algorithm. In this work, the proposed method comprises two main steps. In the first phase, we use the Gauss high pass filter on the star's image after performing bispectrum transformation, and we obtain the direction parameter θ of motion blur by applying PCA. To estimate length L, we rotate the direction of linear motion of the star image to horizontal. en, we use Radon transform to obtain the main lobe of the curve in the vertical direction. Finally, we obtained the length L by finding the first minimum of the curve. In the next phase, we reconstruct the star points in the tracking window, and the pending area is reconstructed by using the hyper-Laplacian prior method and blur kernel.

Blur Kernel Estimation
In this section, we present the method for angle estimation. In the frequency domain of the image, the direct estimation of the angle results in large number of errors due to additive noise. is significantly affects the process of parameter estimation. Generally, the noise added in the star's image is Gaussian noise with zero means, and its higher-order cumulant is zero. erefore, we accomplish the angle estimation in the bispectrum domain [14] of the image. We first perform the bispectrum transformation of the acquired image. is is presented in Figure 4, where Figure 4(b) represents the center area of the bispectrum, and the peak of gray value in the center of the image represents the noise accumulation. It is noticeable that the higher the noise level, the higher is the peak. We use the Gauss high pass filter on the star's image after performing bispectrum transformation. As presented in Figure 4(c), the direction of the ellipse represents the direction of motion. In order to ensure the accuracy and speed of parameter estimation, we obtain the central ellipse region with directional characteristics after performing the threshold segmentation of Figure 4(c). Finally, we obtain the direction parameter θ of motion blur by applying PCA on the central bright region. e algorithm of angle calculation using PCA is as follows Algorithm 1.
After applying threshold segmentation in bispectrum domain, the region pixel is set presented as e total number of pixels is N. Now, P � [p 1 , p 2 . . . n] T , and P is the matrix of dimensions n × 2. Each row in this matrix represents the coordinate value of each pixel. Afterward, we calculate the covariance matrix according to C � (P − μ) T (P − μ), where µ indicates the mean of P. Finally, the eigenvalue decomposition of matrix C is carried out to obtain the eigenvector v max � (w 1 , w 2 ) and the corresponding highest eigenvalue λ max . Now, the motion angle is calculated as

Length Estimation.
e Fourier transform G (u, v) of the star's image with linear motion blur, i.e., g (x, y), seems like a diffraction fringe. Generally, according to the properties of the Fourier transform, the direction of the stripes is perpendicular to the direction of motion. In addition, the spacing of the stripes is dependent on the length of the motion. In order to perform the quantitative and uniform analysis of the relationship between stripe spacing and length of the motion, we rotate the direction of linear motion of star image from vertical to horizontal. e rotated image is expressed as g r (x, y). e Fourier transform of this image is G r (u, v).
In addition, the relative intensity of G r (u, v) along the vertical direction is where N represents the height of the image and L represents the motion length of the motion blur. H(ξ) is an ideal sinc. erefore, there is a relation between L and the width of the main lobe of H(ξ): where d denotes the distance between two adjacent dark stripes, which is also equivalent to half of the main lobe. erefore, the process of obtaining d from G (u, v) is crucial. e Radon transform is an effective tool to perform this computation. As presented on the left side of Figure 5, the energy of main lobe of G r (u, v) is significantly large as compared to the energy of the side lobes. us, we use Radon transform to obtain the main lobe of H(ξ) in the vertical direction. e Radon transform calculates the projection of the image matrix in a specific direction. e Radon transform in the continuous domain is mathematically defined as As discussed, we obtain an estimate of H(ξ) and R(ξ, 0). e right side of Figure 5 exhibits the result of estimation (horizontal Radon transform). For convenience, we denote the estimation of H(ξ) as H e (ξ).
In real-time application, the curve obtained by Radon transform H e (ξ) is not an ideal sinc, due to the addition of noise. e light and dark fringes, i.e., the peaks and valleys of the sinc curve gradually weaken with the increase in the noise level. is is presented in Figure 6. In this work, we propose a robust method to compute the first minimum of the sinc curve. We successfully estimate the length, even under the low SNR and small blur length. e pseudocode of the proposed method is presented in Algorithm 2. Input: the star image g(x, y) (1) Compute the bispectrum transformation of g(x, y), G 2 (u′, v′) (2) Apply Gauss high pass filtering on G 2 (u′, v′), and obtain G H (2) (3) Perform threshold segmentation on G H (2) , and obtain the position matrix of bright region pixels P (4) Compute the mean of P, µ (5) Compute the covariance matrix C ←(P − μ) T (P − μ) (6) Compute the maximum eigenvalue λ max , and its eigenvector v max (7) θ ← arctan(|w 2 /w 1 |) Output: θ Mathematical Problems in Engineering affected by noise. e proposed method estimates whether the d obtained in Step 5 is abnormal. If the estimated parameter is abnormal, we adjust the weight � 1.5 and recalculate the parameter. e curve presented in Figure 6(c) is obtained when the blur length L � 5 and the noise is present. e curve presented in Figure 6(d) presents the curve that is processed without adjusting the weight. e curve presented in Figure 6(e) is processed after adjusting the weight.

Area Filtering.
e global recovery is unable to meet the real-time requirements of the star sensor. erefore, in this work, we propose the area filtering method. We obtain the center position (x, y) of the star point in the tracking window from the previous frame and the motion length L. e area, i.e., (2 L) * (2 L), is selected as the pending area with (x, y) as the center. e pending area is reconstructed by using the hyper-Laplacian prior method and the centroid of star point is extracted.
Equation (9) is the cost function of the hyper-Laplace restoration method [19], where x is the restored image, k is the input blur kernel, and y is the image degraded due to blur and noise. Assume that y is obtained by adding noise after convolving x and k. Given y and x, x is obtained by minimizing the cost function. λ is the parameter that controls the intensity of regularization, and |x ⊕ f i | α is the penalty function. When α � 1/2 and α � 2/3, there is an analytical solution.

Experiments and Results
In this section, we present the evaluation of the proposed algorithm by using both the simulated star maps and the real star maps. We divide the simulations into three phases. In the first phase, we verify that the proposed blur kernel estimation algorithm has the capability to quickly estimate the high-precision blur kernel. In the second phase, we verify that the proposed area filtering method has the ability to remove the star trailing on the basis of estimated motion parameters in a computationally efficient manner. In the last phase, we verify that the star image processed by this algorithm effectively improves the number of stars and the centroid extraction accuracy of these stars.
We use MATLAB to perform the aforementioned experiment. e computer used for computations has an Intel i5−5200 u cpu and 8 GB memory.

Simulated Star Point.
We simulate the star point and convolve it with the blur kernel to obtain the blur star point after motion blur. In Figure 7, we present the two-dimensional and three-dimensional images of the original and blur stars. Note that the blur length L � 20 pixels.

Blur Kernel Estimation.
e traditional method performs Radon transformation on all angles of the frequency spectrum and seeks the maximum value to obtain the angle parameters of the blur kernel. Similarly, the cepstrum method performs Radon transformation on all angles of the cepstrum to obtain angle parameters. e method proposed in this paper does not need to traverse all angles for Radon transformation but uses PCA to directly calculate angles and can effectively suppress zero-mean Gaussian noise in the dual spectrum domain. ere are different methods which are used to estimate the angle. e results are presented in Figures 8(a) and 8(b). e fixed motion blur length is 20 pixels, and the angle of motion blur is changed from 0°to 90°. We compare the error in the angle estimation of the proposed method with the traditional Radon method, cepstrum + Radon method, and the method proposed by Jiang [30]. e results make it evident that the error for the proposed method from 0°to 90°is smaller as compared to the other methods. In the noiseless case, the maximum error of angle estimation is less than 0.5°. In the noisy case, the maximum error is less than 0.8°. It is noticeable that the estimates of the other three methods pose a lot of fluctuations. is is because the Radon transform is easily affected by noise and image resolution. Although Jiang [30] incorporated preprocessing to reduce the impact of noise, the authors use Radon transform to compute the angle, and thus, its accuracy is not high. On the contrary, the proposed method not only uses the bispectrum to suppress the impact of noise but also uses the PCA to estimate the angle. is greatly improves the accuracy of angle estimation.
For length estimation, the traditional method and the cepstrum method perform Radon transformation according to the blur angle direction to obtain the sinc curve and find the distance 2 d between the first minimum values to obtain the blur length L. In order to suppress the degeneration of noise on the curve, we propose an adjustable weighting method, which converts the search for the distance of the first minimum value to the search for the minimum value, thereby robustly estimating the accurate blur length.
In order to estimate the length, we use the same two groups of tests. When the motion blur angle is 20°, the motion blur length of the simulated star map is transformed from 0 to 70 pixels. We use four methods to estimate the blur length L and calculate the error. Figures 9(a) and 9(b) present the test results for two cases, i.e., with and without noise, respectively. Note that as the blur length increases, the errors of the four methods depict an increasing trend. e performance of the method presented by Jiang [30] in the noisy environment is not as efficient as it is in the noiseless experiment. is may be caused by the excessive preprocessing before parameter estimation. e error of the proposed method is smaller than that of other methods. In addition, the proposed method has a better antinoise performance. Under the high dynamic conditions, the star sensor algorithm must be computationally efficient. erefore, we compare the computational cost of the four methods. We test each of the four methods 50 times using the same star maps, as discussed in the previous subsection, and compute the average time consumed by each method. As presented in Table 1, the method proposed by Jiang [30] is the most computationally expensive algorithm. is is mainly because of the additional preprocessing. e proposed method has the lowest computational cost. e average time is 0.0514 s, which is only one tenth of the traditional method and cepstrum method.

Area Filtering Method.
In this section, we evaluate the star point reconstruction of the proposed algorithm. e method in this paper is based on the hyper-Laplace recovery method (parameter settings: λ � 3000, α � 2/3). e results reveal that the proposed algorithm performs efficient star tracking under the condition of low SNR. We use the proposed method and other methods presented in the literature against the same star   Table 2. Figure 10 shows the degraded star map without noise (when the blur length L � 20, and angle � 20°), using the traditional method (traditional spectrum method of blur kernel estimation and Wiener filter to restore the image), Jiang method [30], and the proposed method to estimate the blur kernel and reconstruct the star points from a single frame star map. It is evident from Figure 10 that the proposed method has the best performance. e method proposed by Jiang [30] is the second best performing algorithm. Figure 11 shows the degraded star map after adding noise (SNR � −3.2 dB, PSNR � 119.2794, RMSE � 0.0166, blur length L � 20, and angle � 20°). Figure 11 shows that the Wiener filter has a poor ability to suppress noise. Similarly, the Jiang method is affected by the ringing phenomenon, and the ringing effect increases with the increase in noise level. On the contrary, the method proposed in this work effectively suppresses the noise and reconstructs the star points successfully. Table 2 presents the image indexes after restoration of Figure 11. Figure 12 presents the real star map, and Table 3 presents the parameters of the real star sensor. Note that the Wiener filter fails to reconstruct the star points. e star points reconstructed by the method proposed by Jiang [30] suffer from the ringing phenomenon, which looks like two stars. Moreover, the number of iterations after which the algorithm should be stopped is selected manually. On the contrary, the method proposed in this work automatically stops and reconstructs the star points successfully and effectively.

Performance of Time.
We analyze the computational time of three blind restoration methods that are traditional method, Jiang method, and the proposed method from blur kernel estimation to star reconstruction. e details of the time consumed by different algorithms are presented in Table 4. Although the method proposed by Jiang accelerates the RL algorithm and the recovery time is 0.0861 s, the excessive preprocessing of blur kernel estimation takes 1.6517 s. e total time consumed in the execution of the proposed algorithm is only 0.1138 s. As compared with the traditional method, the processing speed of the proposed method is 5 times greater. Similarly, the proposed method is 15 times faster as compared with the Jiang method.

Star Centroid Extraction.
In order to verify the effectiveness of the proposed algorithm for improving the centroid accuracy of star points, we perform the centroid extraction of star points before and after reconstruction. e results are presented in Table 5. It is notable that the centroid accuracy of the processed star points is significantly improved. e extraction for stars 2, 3, 8, 9, and 10 failed because the gray level of the star region is too low before processing. So, it is unable to reach the local threshold. e processed star points are gathered by energy; thus, they are successfully segmented and extracted by the local threshold. In short, the algorithm proposed in this work increases the number of star point extraction and the accuracy of star centroid under high dynamic conditions.

Conclusions
In this work, we propose a fast blur kernel estimation method and an area filtering method for star sensors. As compared with other methods, the proposed method has higher accuracy, is robust to noise, and is computationally inexpensive. We use the proposed area filtering method to reconstruct the star points by using the estimated blur kernel parameters. e energy of the reconstructed star points is concentrated. is provides an estimate of the original energy distribution of the star points. us, the trailing of the star points is efficiently removed. As compared with the traditional methods presented in the literature, the computational efficiency of the complete algorithm is 15 times superior to that of the accelerated iterative method proposed by Jiang. e number of star points and the accuracy of the centroid are also improved. In future work, in order to further verify and expand the practicability of the proposed algorithm, we will continue to improve the overall real-time performance and accuracy. In addition, we will also perform the application test of the algorithm on the star sensor combined with hardware.

Data Availability
No data were used to support this study..

Conflicts of Interest
e authors declare that they have no conflicts of interest.