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

Authors: Yaxian Hou; Rujin Zhao; Yuebo Ma; Longdong He; Zifa Zhu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The 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. The motion 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. This algorithm quickly reconstructs the star points in the tracking windows and effectively removes the star point tailing in real time. The computational efficiency of the proposed algorithm is 5 times superior to the traditional method and 15 times superior to the existing accelerated iterative method. The 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.


---

## Body

## 1. Introduction

Star sensor device forms the basis of celestial navigation. The three-axis attitude information of spacecraft is calculated by using the observed star’s position [1]. The star sensor has three working modes: whole sky mode, local sky mode, and tracking mode. The tracking mode is the main mode of star sensor’s working.The star sensor performs poorly under high dynamic conditions. This 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 Figure1, the motion blur disperses the energy of stars in the image, consequently, leading to the reduction in the signal-to-noise ratio (SNR). This 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. The aforementioned factors determine if the star sensor is able to perform normally under high dynamic conditions. Therefore, 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.Figure 1
The process of star spot tailing.There are two major techniques that are used to enhance the accuracy of star centroid extraction, namely, hardware enhancement and software restoration. The 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. This 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. Thus, we need additional software processing for achieving our goal.The software restoration method mainly removes the tail of the star spot for improving the accuracy of star’s centroid extraction. The 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–6], Hough transform, and other line detection techniques to estimate the motion blur angle in the frequency domain. The 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. The technique of Radon transform traverses all the angles on the spectrum chart to find the maximum value. This maximum value is the estimated angle, as presented in Figure 2. The Radon transform is not only affected easily by image resolution and noise, but it also has poor robustness and higher computation cost.Figure 2
The Radon transform of all angles.In recent years, deep learning [13–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]. There are also many research studies on restoration methods based on deep learning. The literature [20] combines the advantages of energy functional optimization and deep learning inverse convolution models and uses half-duadratic splitting (HQS) to decompose the energy functional and applies CNN noise prior deals with the subproblem of denoising. The 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. This 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. This 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–26]. In [27], the authors propose a star centroid motion model according to the angular velocity obtained from the gyroscope. This 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. Thus, 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. The 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. This is because this method requires the complete image of the star to be processed, while the star sensor is in the tracking mode. Thus, the algorithm is unable to meet the real-time tracking requirements of the star sensor.The process of reconstructing the image is to use inverse convolution to reverse the image degradation process. The energy functional inverse convolution model depends on the image prior. The difference between the different methods is roughly that the norm used in the model is different. The 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. The 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. Therefore, 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. The 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. Therefore, 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–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. The 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].The rest of the manuscript is organized as follows: in Section2, 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.


## 2. Mechanism Analysis

The 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 distributionE (x, y).(1)Ex,y=E02πσ2exp−x−x02+y−y022σ2,where σ represents the variance of energy distribution (σ is dependent on the parameters of star sensor) and E0 represents the total energy of the stars. x0,y0 presents the center point. The star image f (x, y) is represented by(2)fx,y=C⋅Ex,y,C=const.When the star sensor operates in a high dynamic environment, there is a relative motion between the star and the sensor. The 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]. The energy distribution E (x, y) of the star’s image disperses on the sensor, which is in a relative motion to the star. This phenomenon leads to the formation of “tail” of the star spot. In addition, the center of the star spot is also displaced. The 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(3)Ex,y=E02πσ2L∫0Lexp−Δx2+Δy22σ2dt,where ∆x = x − x0 − rcosα and ∆y = y − y0–rsinα. α, r, and L denote the motion angle, displacement, and star point motion length in the exposure time δt, respectively.The blurred image is represented asg(x, y) (g(x, y) = C·Eb). 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:(4)gx,y=fx,y⊗hx,y+nx,y,where h (x, y) represents the motion blur kernel. This kernel contains the motion information and interacts with the noise in the image to get the degraded image.


## 3. A Star Tailing Removal Method Based on Fast Blur Kernel Estimations

Figure3 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. Then, 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.Figure 3
The block diagram for the proposed algorithm.


### 3.1. Blur Kernel Estimation

#### 3.1.1. Angle 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. This 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. Therefore, we accomplish the angle estimation in the bispectrum domain [14] of the image. We first perform the bispectrum transformation of the acquired image. This 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.Figure 4
Process of angle estimation: (a) the three-dimensional display of the bispectrum of the blurred image, (b) the bispectrum, and (c) the bispectrum after Gauss high pass filtering.
(a)(b)(c)The algorithm of angle calculation using PCA is as follows Algorithm1.Algorithm 1: Angle estimation based on PCA.
Input: the star image g(x, y)(1)
Compute the bispectrum transformation ofg(x, y),G2u′,v′(2)
Apply Gauss high pass filtering onG2u′,v′, and obtain GH2(3)
Perform threshold segmentation onGH2, and obtain the position matrix of bright region pixels P(4)
Compute the mean ofP,µ(5)
Compute the covariance matrixC ←P−μTP−μ(6)
Compute the maximum eigenvalueλmax, and its eigenvector vmax(7)
θ←arctan|w2/w1|Output:θAfter applying threshold segmentation in bispectrum domain, the region pixel is set presented asPk=xk,ykT. The total number of pixels is N. Now, P=p1,p2…nT, 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−μTP−μ, where µ indicates the mean of P. Finally, the eigenvalue decomposition of matrix C is carried out to obtain the eigenvector vmax=w1,w2 and the corresponding highest eigenvalue λmax. Now, the motion angle is calculated as(5)θ=arctanw2w1.


#### 3.1.2. Length Estimation

The Fourier transformG (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. The rotated image is expressed as grx,y. The Fourier transform of this image is Gru,v.In addition, the relative intensity ofGru,v along the vertical direction is(6)Hξ=SincLξN,ξ=0,1,2…..N−1,where N represents the height of the image and L represents the motion length of the motion blur. Hξ is an ideal sinc. Therefore, there is a relation between L and the width of the main lobe of Hξ:(7)L=Nd,where d denotes the distance between two adjacent dark stripes, which is also equivalent to half of the main lobe. Therefore, the process of obtaining d from G (u, v) is crucial.The Radon transform is an effective tool to perform this computation. As presented on the left side of Figure5, the energy of main lobe of Gru,v is significantly large as compared to the energy of the side lobes. Thus, we use Radon transform to obtain the main lobe of Hξ in the vertical direction. The Radon transform calculates the projection of the image matrix in a specific direction. The Radon transform in the continuous domain is mathematically defined as(8)Rρ,θ=∫−∞+∞Gru,vδρ−ucosθ−vsinθdudv.Figure 5
The image ofGru,v (on the right is Heξ obtained by Randon transform).As discussed, we obtain an estimate ofHξ and Rξ,0. The right side of Figure 5 exhibits the result of estimation (horizontal Radon transform). For convenience, we denote the estimation of Hξ as Heξ.In real-time application, the curve obtained by Radon transformHeξ is not an ideal sinc, due to the addition of noise. The light and dark fringes, i.e., the peaks and valleys of the sinc curve gradually weaken with the increase in the noise level. This 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. The pseudocode of the proposed method is presented in Algorithm 2.Figure 6
The estimation ofHeξ with noise: (a) sinc curve (when L = 20), (b) processed curve, (c) sinc curve (when L = 5), (d) processed result without weight adjustment, and (e) processed result after weight adjustment.
(a)(b)(c)(d)(e)Algorithm 2: The proposed algorithm for length estimation.
Input:Hξ(1)
Compute the length ofHξ, N(2)
β←4(3)
Pi←Hξ1(4)
fori = 1 to 2 do(5)
forj = 1 toN-1 do(6)
Δj←Hξj+1−Hξj(7)
ifΔj>0then(8)
Δj′←βΔj(9)
else(10)
Δj′←Δj(11)
end if(12)
Pj+1←Pi+Δj′(13)
end for(14)
d←minP(15)
ifd<5then(16)
β←1.5(17)
else(18)
break(19)
end if(20)
end for(21)
L←N/dOutput: LFigure6(a) represents Heξ, when the blur length L = 20 with noise, and Figure 6(b) represents the curve processed by this method. Since the blur length L is very small, it is easily affected by noise. The 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. The curve presented in Figure 6(c) is obtained when the blur length L = 5 and the noise is present. The curve presented in Figure 6(d) presents the curve that is processed without adjusting the weight. The curve presented in Figure 6(e) is processed after adjusting the weight.


### 3.2. Area Filtering

The global recovery is unable to meet the real-time requirements of the star sensor. Therefore, 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. The area, i.e., (2 L) ∗ (2 L), is selected as the pending area with (x, y) as the center. The 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⊕fiα is the penalty function. When α = 1/2 and α = 2/3, there is an analytical solution.(9)min∑i=1Nλ2x⊕k−yi2+∑j=1Jx⊕fiα.


## 3.1. Blur Kernel Estimation

### 3.1.1. Angle 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. This 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. Therefore, we accomplish the angle estimation in the bispectrum domain [14] of the image. We first perform the bispectrum transformation of the acquired image. This 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.Figure 4
Process of angle estimation: (a) the three-dimensional display of the bispectrum of the blurred image, (b) the bispectrum, and (c) the bispectrum after Gauss high pass filtering.
(a)(b)(c)The algorithm of angle calculation using PCA is as follows Algorithm1.Algorithm 1: Angle estimation based on PCA.
Input: the star image g(x, y)(1)
Compute the bispectrum transformation ofg(x, y),G2u′,v′(2)
Apply Gauss high pass filtering onG2u′,v′, and obtain GH2(3)
Perform threshold segmentation onGH2, and obtain the position matrix of bright region pixels P(4)
Compute the mean ofP,µ(5)
Compute the covariance matrixC ←P−μTP−μ(6)
Compute the maximum eigenvalueλmax, and its eigenvector vmax(7)
θ←arctan|w2/w1|Output:θAfter applying threshold segmentation in bispectrum domain, the region pixel is set presented asPk=xk,ykT. The total number of pixels is N. Now, P=p1,p2…nT, 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−μTP−μ, where µ indicates the mean of P. Finally, the eigenvalue decomposition of matrix C is carried out to obtain the eigenvector vmax=w1,w2 and the corresponding highest eigenvalue λmax. Now, the motion angle is calculated as(5)θ=arctanw2w1.


### 3.1.2. Length Estimation

The Fourier transformG (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. The rotated image is expressed as grx,y. The Fourier transform of this image is Gru,v.In addition, the relative intensity ofGru,v along the vertical direction is(6)Hξ=SincLξN,ξ=0,1,2…..N−1,where N represents the height of the image and L represents the motion length of the motion blur. Hξ is an ideal sinc. Therefore, there is a relation between L and the width of the main lobe of Hξ:(7)L=Nd,where d denotes the distance between two adjacent dark stripes, which is also equivalent to half of the main lobe. Therefore, the process of obtaining d from G (u, v) is crucial.The Radon transform is an effective tool to perform this computation. As presented on the left side of Figure5, the energy of main lobe of Gru,v is significantly large as compared to the energy of the side lobes. Thus, we use Radon transform to obtain the main lobe of Hξ in the vertical direction. The Radon transform calculates the projection of the image matrix in a specific direction. The Radon transform in the continuous domain is mathematically defined as(8)Rρ,θ=∫−∞+∞Gru,vδρ−ucosθ−vsinθdudv.Figure 5
The image ofGru,v (on the right is Heξ obtained by Randon transform).As discussed, we obtain an estimate ofHξ and Rξ,0. The right side of Figure 5 exhibits the result of estimation (horizontal Radon transform). For convenience, we denote the estimation of Hξ as Heξ.In real-time application, the curve obtained by Radon transformHeξ is not an ideal sinc, due to the addition of noise. The light and dark fringes, i.e., the peaks and valleys of the sinc curve gradually weaken with the increase in the noise level. This 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. The pseudocode of the proposed method is presented in Algorithm 2.Figure 6
The estimation ofHeξ with noise: (a) sinc curve (when L = 20), (b) processed curve, (c) sinc curve (when L = 5), (d) processed result without weight adjustment, and (e) processed result after weight adjustment.
(a)(b)(c)(d)(e)Algorithm 2: The proposed algorithm for length estimation.
Input:Hξ(1)
Compute the length ofHξ, N(2)
β←4(3)
Pi←Hξ1(4)
fori = 1 to 2 do(5)
forj = 1 toN-1 do(6)
Δj←Hξj+1−Hξj(7)
ifΔj>0then(8)
Δj′←βΔj(9)
else(10)
Δj′←Δj(11)
end if(12)
Pj+1←Pi+Δj′(13)
end for(14)
d←minP(15)
ifd<5then(16)
β←1.5(17)
else(18)
break(19)
end if(20)
end for(21)
L←N/dOutput: LFigure6(a) represents Heξ, when the blur length L = 20 with noise, and Figure 6(b) represents the curve processed by this method. Since the blur length L is very small, it is easily affected by noise. The 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. The curve presented in Figure 6(c) is obtained when the blur length L = 5 and the noise is present. The curve presented in Figure 6(d) presents the curve that is processed without adjusting the weight. The curve presented in Figure 6(e) is processed after adjusting the weight.


## 3.1.1. Angle 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. This 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. Therefore, we accomplish the angle estimation in the bispectrum domain [14] of the image. We first perform the bispectrum transformation of the acquired image. This 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.Figure 4
Process of angle estimation: (a) the three-dimensional display of the bispectrum of the blurred image, (b) the bispectrum, and (c) the bispectrum after Gauss high pass filtering.
(a)(b)(c)The algorithm of angle calculation using PCA is as follows Algorithm1.Algorithm 1: Angle estimation based on PCA.
Input: the star image g(x, y)(1)
Compute the bispectrum transformation ofg(x, y),G2u′,v′(2)
Apply Gauss high pass filtering onG2u′,v′, and obtain GH2(3)
Perform threshold segmentation onGH2, and obtain the position matrix of bright region pixels P(4)
Compute the mean ofP,µ(5)
Compute the covariance matrixC ←P−μTP−μ(6)
Compute the maximum eigenvalueλmax, and its eigenvector vmax(7)
θ←arctan|w2/w1|Output:θAfter applying threshold segmentation in bispectrum domain, the region pixel is set presented asPk=xk,ykT. The total number of pixels is N. Now, P=p1,p2…nT, 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−μTP−μ, where µ indicates the mean of P. Finally, the eigenvalue decomposition of matrix C is carried out to obtain the eigenvector vmax=w1,w2 and the corresponding highest eigenvalue λmax. Now, the motion angle is calculated as(5)θ=arctanw2w1.


## 3.1.2. Length Estimation

The Fourier transformG (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. The rotated image is expressed as grx,y. The Fourier transform of this image is Gru,v.In addition, the relative intensity ofGru,v along the vertical direction is(6)Hξ=SincLξN,ξ=0,1,2…..N−1,where N represents the height of the image and L represents the motion length of the motion blur. Hξ is an ideal sinc. Therefore, there is a relation between L and the width of the main lobe of Hξ:(7)L=Nd,where d denotes the distance between two adjacent dark stripes, which is also equivalent to half of the main lobe. Therefore, the process of obtaining d from G (u, v) is crucial.The Radon transform is an effective tool to perform this computation. As presented on the left side of Figure5, the energy of main lobe of Gru,v is significantly large as compared to the energy of the side lobes. Thus, we use Radon transform to obtain the main lobe of Hξ in the vertical direction. The Radon transform calculates the projection of the image matrix in a specific direction. The Radon transform in the continuous domain is mathematically defined as(8)Rρ,θ=∫−∞+∞Gru,vδρ−ucosθ−vsinθdudv.Figure 5
The image ofGru,v (on the right is Heξ obtained by Randon transform).As discussed, we obtain an estimate ofHξ and Rξ,0. The right side of Figure 5 exhibits the result of estimation (horizontal Radon transform). For convenience, we denote the estimation of Hξ as Heξ.In real-time application, the curve obtained by Radon transformHeξ is not an ideal sinc, due to the addition of noise. The light and dark fringes, i.e., the peaks and valleys of the sinc curve gradually weaken with the increase in the noise level. This 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. The pseudocode of the proposed method is presented in Algorithm 2.Figure 6
The estimation ofHeξ with noise: (a) sinc curve (when L = 20), (b) processed curve, (c) sinc curve (when L = 5), (d) processed result without weight adjustment, and (e) processed result after weight adjustment.
(a)(b)(c)(d)(e)Algorithm 2: The proposed algorithm for length estimation.
Input:Hξ(1)
Compute the length ofHξ, N(2)
β←4(3)
Pi←Hξ1(4)
fori = 1 to 2 do(5)
forj = 1 toN-1 do(6)
Δj←Hξj+1−Hξj(7)
ifΔj>0then(8)
Δj′←βΔj(9)
else(10)
Δj′←Δj(11)
end if(12)
Pj+1←Pi+Δj′(13)
end for(14)
d←minP(15)
ifd<5then(16)
β←1.5(17)
else(18)
break(19)
end if(20)
end for(21)
L←N/dOutput: LFigure6(a) represents Heξ, when the blur length L = 20 with noise, and Figure 6(b) represents the curve processed by this method. Since the blur length L is very small, it is easily affected by noise. The 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. The curve presented in Figure 6(c) is obtained when the blur length L = 5 and the noise is present. The curve presented in Figure 6(d) presents the curve that is processed without adjusting the weight. The curve presented in Figure 6(e) is processed after adjusting the weight.


## 3.2. Area Filtering

The global recovery is unable to meet the real-time requirements of the star sensor. Therefore, 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. The area, i.e., (2 L) ∗ (2 L), is selected as the pending area with (x, y) as the center. The 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⊕fiα is the penalty function. When α = 1/2 and α = 2/3, there is an analytical solution.(9)min∑i=1Nλ2x⊕k−yi2+∑j=1Jx⊕fiα.


## 4. 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. The computer used for computations has an Intel i5−5200 u cpu and 8 GB memory.


### 4.1. 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 Figure7, we present the two-dimensional and three-dimensional images of the original and blur stars. Note that the blur length L = 20 pixels.Figure 7
The simulated star point: (a) simulated star image and (b) simulated trailing star image.
(a)(b)


### 4.2. Blur Kernel Estimation

The 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. The 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.There are different methods which are used to estimate the angle. The results are presented in Figures8(a) and 8(b). The 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]. The 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. This 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. This greatly improves the accuracy of angle estimation.Figure 8
The estimation error of different angles: (a) noiseless case and (b) noisy case (SNR = 0.4979 dB).
(a)(b)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 2d 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 lengthL 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. The performance of the method presented by Jiang [30] in the noisy environment is not as efficient as it is in the noiseless experiment. This may be caused by the excessive preprocessing before parameter estimation. The error of the proposed method is smaller than that of other methods. In addition, the proposed method has a better antinoise performance.Figure 9
The estimation errors of different lengths: (a) noiseless case and (b) noisy case (SNR = 0.4979 dB).
(a)(b)Under the high dynamic conditions, the star sensor algorithm must be computationally efficient. Therefore, 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 Table1, the method proposed by Jiang [30] is the most computationally expensive algorithm. This is mainly because of the additional preprocessing. The proposed method has the lowest computational cost. The average time is 0.0514 s, which is only one tenth of the traditional method and cepstrum method.Table 1
The comparison of computational complexity of different methods.
Time average (s)Proposed method0.0514Traditional Radon transform method0.4350Cepstrum + Radon transform1.6517


### 4.3. Area Filtering Method

In this section, we evaluate the star point reconstruction of the proposed algorithm. The method in this paper is based on the hyper-Laplace recovery method (parameter settings:λ = 3000, α = 2/3). The 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 map. The results are presented in Figures 10–12 and 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. The 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 10
The noise-free reconstruction of simulated star points: (a) degenerated star points, (b) traditional method, (c) method proposed by Jiang [30], and (d) the proposed method.
(a)(b)(c)(d)Figure 11
The reconstruction of simulated star with noise (SNR = −3.2 dB, PSNR = 119.2794, RMSE = 0.0166,L = 20, and angle = 20°): (a) degenerated star points, (b) traditional method, (c) method proposed by Jiang [30], and (d) the proposed method.
(a)(b)(c)(d)Figure 12
The reconstruction of real stars: (a) degenerated star points, (b) traditional method, (c) method proposed by Jiang [30], and (d) the proposed method.
(a)(b)(c)(d)Table 2
The results after processing by different methods.
Traditional method [26]Jiang method [30]Proposed methodBlur kernelAngle232120.2141Length18.692319.015120.0784SNR/dB−1.06971.54675.2078PSNR/dB117.6414122.8742127.7242RMSE0.01830.01350.0089Figure12 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. The 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.Table 3
The parameters of real star sensor.
Star sensor parametersFOV (°)20∗20Time of exposure (ms)118Frame rate (Hz)10Imager size (pixel)1024∗1024


### 4.4. 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. The details of the time consumed by different algorithms are presented in Table4. 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. The 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.Table 4
The details of time consumption with different methods.
Blur kernel estimationReconstructTotal timeTraditional method0.43500.05450.4895Jiang method1.65170.08611.7334Proposed method0.05140.06240.1138


### 4.5. 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.The results are presented in Table5. It is notable that the centroid accuracy of the processed star points is significantly improved. The 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. The 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.Table 5
The error of center of mass of stars before and after processing (absolute value).
NoBeforeAfterΔx/pixelΔy/pixelΔx/pixelΔy/pixel10.11620.15320.02000.04012Extraction failed0.00580.027623Extraction failed0.03790.055740.21540.10070.01020.066150.27940.18310.06880.161460.29010.31040.29150.121575.78732.07710.42750.09228Extraction failed0.10170.04959Extraction failed0.11040.059210Extraction failed0.47400.1684


## 4.1. 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 Figure7, we present the two-dimensional and three-dimensional images of the original and blur stars. Note that the blur length L = 20 pixels.Figure 7
The simulated star point: (a) simulated star image and (b) simulated trailing star image.
(a)(b)


## 4.2. Blur Kernel Estimation

The 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. The 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.There are different methods which are used to estimate the angle. The results are presented in Figures8(a) and 8(b). The 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]. The 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. This 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. This greatly improves the accuracy of angle estimation.Figure 8
The estimation error of different angles: (a) noiseless case and (b) noisy case (SNR = 0.4979 dB).
(a)(b)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 2d 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 lengthL 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. The performance of the method presented by Jiang [30] in the noisy environment is not as efficient as it is in the noiseless experiment. This may be caused by the excessive preprocessing before parameter estimation. The error of the proposed method is smaller than that of other methods. In addition, the proposed method has a better antinoise performance.Figure 9
The estimation errors of different lengths: (a) noiseless case and (b) noisy case (SNR = 0.4979 dB).
(a)(b)Under the high dynamic conditions, the star sensor algorithm must be computationally efficient. Therefore, 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 Table1, the method proposed by Jiang [30] is the most computationally expensive algorithm. This is mainly because of the additional preprocessing. The proposed method has the lowest computational cost. The average time is 0.0514 s, which is only one tenth of the traditional method and cepstrum method.Table 1
The comparison of computational complexity of different methods.
Time average (s)Proposed method0.0514Traditional Radon transform method0.4350Cepstrum + Radon transform1.6517


## 4.3. Area Filtering Method

In this section, we evaluate the star point reconstruction of the proposed algorithm. The method in this paper is based on the hyper-Laplace recovery method (parameter settings:λ = 3000, α = 2/3). The 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 map. The results are presented in Figures 10–12 and 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. The 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 10
The noise-free reconstruction of simulated star points: (a) degenerated star points, (b) traditional method, (c) method proposed by Jiang [30], and (d) the proposed method.
(a)(b)(c)(d)Figure 11
The reconstruction of simulated star with noise (SNR = −3.2 dB, PSNR = 119.2794, RMSE = 0.0166,L = 20, and angle = 20°): (a) degenerated star points, (b) traditional method, (c) method proposed by Jiang [30], and (d) the proposed method.
(a)(b)(c)(d)Figure 12
The reconstruction of real stars: (a) degenerated star points, (b) traditional method, (c) method proposed by Jiang [30], and (d) the proposed method.
(a)(b)(c)(d)Table 2
The results after processing by different methods.
Traditional method [26]Jiang method [30]Proposed methodBlur kernelAngle232120.2141Length18.692319.015120.0784SNR/dB−1.06971.54675.2078PSNR/dB117.6414122.8742127.7242RMSE0.01830.01350.0089Figure12 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. The 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.Table 3
The parameters of real star sensor.
Star sensor parametersFOV (°)20∗20Time of exposure (ms)118Frame rate (Hz)10Imager size (pixel)1024∗1024


## 4.4. 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. The details of the time consumed by different algorithms are presented in Table4. 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. The 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.Table 4
The details of time consumption with different methods.
Blur kernel estimationReconstructTotal timeTraditional method0.43500.05450.4895Jiang method1.65170.08611.7334Proposed method0.05140.06240.1138


## 4.5. 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.The results are presented in Table5. It is notable that the centroid accuracy of the processed star points is significantly improved. The 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. The 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.Table 5
The error of center of mass of stars before and after processing (absolute value).
NoBeforeAfterΔx/pixelΔy/pixelΔx/pixelΔy/pixel10.11620.15320.02000.04012Extraction failed0.00580.027623Extraction failed0.03790.055740.21540.10070.01020.066150.27940.18310.06880.161460.29010.31040.29150.121575.78732.07710.42750.09228Extraction failed0.10170.04959Extraction failed0.11040.059210Extraction failed0.47400.1684


## 5. 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. The energy of the reconstructed star points is concentrated. This provides an estimate of the original energy distribution of the star points. Thus, 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. The 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.



---

*Source: 8819277-2021-02-15.xml*