Restoration of motion-blurred star images with elliptical star streaks

The restoration of motion-blurred star images under high dynamic conditions is important for the high-precision attitude measurement of star sensors. Through motion modelling analysis, it is found that the streak of the imaged star point (star streak) is an elliptical arc. However, existing star image restoration methods are only suitable for the case where the star streak is a straight line. For this reason, a star image restoration algorithm for elliptical star streaks is proposed in this paper. First, the elliptical star streak is transformed into a circular star streak by projective transformation. Then, the circular star streak is transformed into a straight star streak by polar coordinate transformation. Finally, the restored original star image is obtained by restoration methods for straight star streaks and coordinate inverse transformation. At the same time, the algorithm is further optimized by subdividing the polar coordinates. The experiment shows that the proposed algorithm is effective and the restoration accuracy is at the same level as that of existing star image restoration methods for straight star streaks.


Introduction
A star sensor is an instrument that performs high-precision attitude measurements through observations of stars. It first takes images of stars and then performs attitude calculations of the spacecraft based on the imaging information. When the spacecraft is moving at a high speed, the imaged stars will produce a streak which is a kind of motion blur. The star streak will be different under different motion situations. At the same time, the starlight energy will be dispersed to multiple pixels, resulting in a lower image signal-to-noise ratio and a lower star extraction accuracy. This further leads to a reduction in * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the accuracy of the attitude measurement. Therefore, how to eliminate the blur of imaged stars under high dynamic conditions is of great importance for the attitude measurement of the star sensor.
To solve the above motion blur problem, one kind of approach is to mitigate the degradation of imaging quality by improving the hardware, such as Servo Tracking Platform Technology [1], which avoids motion blur of the image by controlling the platform to compensate the movement of the star sensor. Or, image blur can be reduced by sensitivity improved image sensor, such as electron multiplying CCD [2] and intensified CCD [3]. The improvement of sensitivity can reduce the exposure time and thus reduce the image blurring. In addition, by increasing the number of field of view (FOV) of the star sensor, the sensitivity requirement for a single FOV is reduced, so that the exposure time of each FOV can be reduced to suppress image blur, such as the star sensor with multiple FOV designed in [4]. These hardware based methods are relatively difficult and expensive.
The other kind of more commonly used approach is to eliminate motion blur by software. The software-based approach can be sub-divided into two types. One is non-blind restoration. Motion-blurred star images are restored with the help of other instruments which can measure the angular velocity of the spacecraft. Xiaojuan and Xinlong [5] and Fei et al [6] restore motion-blurred star images with the help of strap-down inertial navigation system or MEMS (micro electro mechanical system) gyroscopes. In order to improve the accuracy of spacecraft angular velocity estimation, Chen et al [7] further adds extended Kalman filter to the above systems. Sun et al [8] combines the star sensor with inertial navigation system based on the work of Chen et al [7]. Ma [9] proposes a multiseeded star region growth for star image restoration aided by MEMS gyroscope. All of these methods realize image restoration based on blur kernel, which is generated according to the angular velocity of spacecraft measured by other instruments. In addition, He et al proposed a star image restoration method based on multi-frame superposition. In this method, adjacent star images are corrected by a motion recursive model. Then, they are superimposed together to realize the restoration. This method focuses on the noise removal and quality of restored star image, and it also requires the aid of other measuring instruments.
The other type of software-based approach is blind restoration. The blind restoration is based only on the image taken by the star sensor. Therefore, it is not necessary to know any motion information of the star sensor before the restoration. For the blind restoration of straight star streaks, existing methods first obtain the direction of the motion and the length of the streak. Radon transform [10], Hough transform [11], constructing differential operator [12], curve fitting [13], spectral method [14], bispectrum [15], etc are commonly used methods to obtain the above information. Then the fuzzy kernel function is generated according to that. Finally, image restoration algorithms like L-R filtering [16], Wiener filtering [17], inverse filtering [18], etc are used to restore the motionblurred star image of straight star streaks. In addition, there are a series of optimized and improved algorithms for detecting the angle or length of straight star streaks and restoring images [19][20][21][22][23][24][25][26][27][28][29][30][31]. Representative work is as follows. In [22,23], Huang, Jie 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, this method requires data of an entire image, thus it cannot be applied to star sensor working in tracking mode. In [28], Hou et al use principal component analysis to estimate the angle of blur kernel. In addition, an adjustable weighting method is proposed to estimate the length of blur kernel. So, the method is able to quickly estimate the high-precision blur kernel based on single degraded image. It solves the problem that Randon transform is vulnerable to noise and has poor robustness, so the method has better effect for images with low signal-to-noise ratio.
In recent years, deep learning and neural network have also been applied in the field of image restoration and have achieved excellent results. Representatively, in [30], Agarwal et al proposed a depth expanded RL algorithm by combining classical RL filtering with the radial basis function neural network. It overcomes the disadvantage that the traditional RL algorithm needs to manually set the number of iterations. In [31], Chen et al proposed a method based on sparse representation, hyper-Laplacian priors, and ensemble neural network, which improves the efficiency and robustness of blur kernel parameter estimation. Image restoration based on deep learning can achieve simultaneous estimation of blur kernel and sharp images, but it requires a large number of training samples and a reasonable neural network model.
Blind restoration algorithms described above for star images are only suitable for straight star streaks. However, when the motion situation of the star sensor becomes more complicated, especially when none of the angular velocities of the three axes is 0, the star streak is not a straight line. None of the above methods can be applied.
Meanwhile, through motion modeling analysis, it is found that the star streak is an elliptical arc. Thus, a blind restoration algorithm for elliptical star streaks is proposed in this paper. It provides a new way for the restoration of motion-blurred star images when the star sensor is under triaxial compound motion.
The rest of this paper is arranged as follows. The mathematical model used in subsequent algorithms are described in section 2. Detailed steps of image restoration algorithm in this paper are described in section 3. The performance of this algorithm is verified by simulation and experiment, and the results are illustrated in section 4. The conclusions are given in section 5. Figure 1 shows the coordinate system O − xyz of the star sensor. The relationship between the star point's observation vector w s and its imaging point P on the image sensor is shown in equation (1). Where (x, y) is the coordinate of the centroid of P and f is the imaging focal length of the star sensor,

Imaging model of a star sensor
When the star sensor works under dynamic conditions, the imaging point P will be blurred. And the star streak of P will be different under different motion situations of the star sensor. Assuming that the change in the vector w s from moment t 0 to t 0 + ∆t is represented by the attitude transformation matrix A t0+∆t t0 , then A t0+∆t t0 can be expressed by the triaxial compound angular velocity of the star sensor as in equation (2).
From figure 1 and equation (3), it can be seen that: When ω z = 0, the star streak is a straight line whose direction is related to that of the compound angular velocity of ω x and ω y .
When ω x = 0 and ω y = 0, that is, when the star sensor only moves around the Z axis, the star streak is a circular arc around the Z axis.
When ω x ̸ = 0, ω y ̸ = 0, and ω z ̸ = 0, the star streak is neither a straight line nor a circular arc, but is a more complex arc.
Existing restoration methods of motion-blurred star images are limited to the case of straight star streaks. In case of complex arc star streaks, new solutions should be found. Firstly, a detailed analysis of the star streak under non-zero triaxial angular velocity is presented below. The analysis provides a theoretical basis for the restoration of motion-blurred star images under composite motion conditions.

Imaging model of a star sensor under compound motion
When the star sensor rotates around the X, Y, Z axes with a triaxial compound angular velocity ⃗ ω, since ∆t ≪ f, The motion model of the star's centroid expressed in equation (3) can be approximated as equation (4): Meanwhile, the direction of ⃗ ω is not the normal vector of the imaging plane as shown in figure 2. It can be seen that the imaging plane intercepts the conical surface at an angle of β as a result of the star streak. Therefore, the star streak is an elliptical arc l [25]. Where the conic surface has the optical center O as its vertex, the direction of the star vector as its generatrix, and the angular velocity vector ⃗ ω as its rotation axis. The angle between the generatrix and the rotation axis is denoted as α, Equation (5) defines the eccentricity of an elliptical star streak. From the definition, it can be seen that α is different for different stars, and β is different for different triaxial compound angular velocities. Therefore, different stars with different compound angular velocities will form elliptical star streaks with different eccentricities. Figure 3 shows the different elliptical star streaks formed by different stars at the same compound angular velocity ⃗ ω.

Restoration of star images with elliptical star streaks
Through the above analysis, it is found that the star streak is an elliptical arc under non-zero triaxial angular velocity. Therefore, this paper proposes a coordinate transformation method based on projective transformation according to the star sensor's imaging model. First, elliptical star streaks with different eccentricities are transferred to circular star streaks with eccentricities of 0. Then, the circular star streaks are transformed to straight star streaks by polar coordinate transformation. After that, these straight star streaks can be restored by the existing image restoration algorithm. At last, the image is returned to its original coordinate after two times of coordinate inverse transformation. Consequently, motion-blurred star image with elliptical star streak is restored.
In the process of polar coordinate transformation, if the quantization units of the rotation angle and the rotation radius are not chosen properly, it will lead to the deformation of the restored star spot. This paper solves the problem by subdividing the quantization units of the rotation angle and selecting the optimal quantization unit of the rotation angle.
The detailed flow chart of the algorithm is shown in figure 4.

Ellipse fitting of the star streak
The projective transformation can only be performed after the expression of the elliptical star streak is obtained through ellipse fitting. The general expression of an ellipse is shown in equation (6): In the process of ellipse fitting, the objective function can be defined as the algebraic distance from the input coordinate point (x i , y i ) to the elliptical star streak, as shown in equation (7): We define the vectors M = [A, B, C, D, E, F] T and X = [ ] T , then equation (7) can be written as equation (8): Thus, the ellipse fitting problem can be reduced to one that can be solved using the least squares method, as shown in equation (9):

Transforming elliptical star streaks to circular star streaks based on projective transformation
From the analysis above, it can be seen that when the star sensor is under triaxial compound motion, the imaging star streak is an elliptical arc obtained by intercepting the conic surface with the imaging plane at a certain angle. The eccentricity of the elliptical star streak is different for different angular velocities of the triaxial compound motion and different positions of the star. As shown in figure 5, the imaging star streak l 1 of star 1 in the imaging plane is an elliptical star streak. If the imaging star streak is projected onto a plane π with ⃗ ω as the normal vector, the projected star streak l ′ 1 will become a circular star streak with an eccentricity of 0. At the same time, streaks of other stars in the projection plane π will also become circular star streaks, such as l ′ 2 ,l ′ 3 in figure 5. Therefore, it is only necessary to determine the projective transformation matrix from the imaging plane to its corresponding plane π. Then the transformation of an elliptical star streak to a circular star streak can be realized.
The transformation matrix from the imaging plane to the projection plane π can be expressed by projective transformation matrix H with nine parameters. The relationship between the points on the imaging plane and the points on the projection plane π can be expressed as follows:  be obtained so that the matrix variance of H can be reduced to eight. Therefore, the matrix can be solved by simply taking the coordinates of four pairs of points on any elliptical star streak in the imaging plane and on its corresponding circular star streak in the projection plane π. In order to improve the accuracy, it is better to select these four groups of corresponding points on a star streak away from the center of the image, as l 3 in figure 5. Meanwhile, for the accuracy of the calculation, four special points on the elliptical star streak in the imaging plane is preferred. As shown in figure (10) can be solved by the following method. By substituting each pair of corresponding points (x i , y i ) and (x ′ i , y ′ i ) (i = 1, 2, 3, 4) into equation (10), we can obtain equations (11) and (12): By expanding equations (11) and (12), we can obtain equations (13) and (14): Equations (13) and (14) can be further expressed in the form of Ah = 0 as shown in equation (15): By combining the above equations obtained from four pairs of corresponding points, we can rewrite the equation in the form of Ah = 0 as shown in equation (16): . .
The sum of square error of equation (16) can be expressed as equation (17): By computing the partial derivative of f for h and then letting its value be 0, we can obtain equation (18): For the eigen-decomposition result of A T A, the exact value of h should be equal to the eigenvector corresponding to the eigenvalue whose value is 0. Considering the noise, h can be estimated as the eigenvector corresponding to the eigenvalue closest to 0. Thus, SVD (singular value decomposition) can be used for decomposing the matrix A: where V is the eigenvector matrix of A T A. The eigenvector corresponding to the eigenvalue closest to 0 is the ninth column of V, as shown in equation (20): Finally, we then obtain the projective transformation matrix H from h. Based on this, the inverse matrix of H is obtained as shown in equation (21). It will be used for subsequent coordinate inverse transformation:

Transforming circular star streaks to straight star streaks based on polar coordinate transformation
After transforming to a circular star streak, a polar coordinate transformation is still needed to transform the circular star streak to a straight one.

Transforming circular star streaks to straight star streaks
based on polar coordinate transformation. The circular star streak in the rectangular coordinate system can be transformed into a straight star streak in the polar coordinate system by polar coordinate transformation [26,27]. The polar coordinate system is defined as a system in which the horizontal coordinate is the angle of rotation θ and the vertical coordinate is  (22). In general, the rotation angle is quantified in the unit of 1 • and the rotation radius is quantified in the unit of 1 pixel, (22) In the polar coordinate system, the length of the straight star streak represents the rotation angle of the circular star streak in the rectangular coordinate system. In a time of ∆t, when different stars rotated and blurred around the Z axis, the rotation angle is same but the length of the circular star streak on the imaging plane is different. Therefore, after the polar coordinate transformation, lengths of the transformed straight star streaks corresponding to all circular star streaks are the same, as shown in figure 6.

Blind restoration of motion-blurred images in the polar
coordinate system. After elliptical star streaks being transformed into straight star streaks, it is possible to restore straight star streak through existing methods such as Lucy-Richardson algorithm.
The Lucy-Richardson algorithm is an image restoration algorithm based on the fuzzy kernel function. The algorithm assumes that the pixels in the image obey the Poisson Distribution. Then the maximum likelihood estimation is applied to iteratively compute the degraded image to obtain a clear image that is closest to the original image. The image restoration model can be expressed as equation (23), where f k (x, y) is the estimation result in the kth iteration of the original clear image, f k+1 (x, y) is the estimation result in the (k + 1) iteration, g(x, y) is the degraded image and h(x, y) is the fuzzy kernel function. The initial value of iterations can be set as f 0 (x, y) = g(x, y). As the number of iterations increases, f k+1 (x, y) will converge gradually, until the original clear image is restored, . (23) To generate the fuzzy kernel function h(x, y), the angle and the length of the motion blur has to be determined first. The angle is 0 for the transformed horizontal straight star streaks, thus only the motion length needs to be acquired. In this paper, the Steger algorithm is used to get the length of the motion blur.
For the motion blur caused by uniform linear motion, the fuzzy kernel function can be expressed as equation (24), where θ is the motion blur angle (the angle between the motion direction and the positive direction of the horizontal axis), L is the motion blur scale (the distance of the pixel movement in the motion direction), Figure 7 shows the restoration results of straight star streaks in the polar coordinate system using the Lucy-Richardson algorithm.

Coordinate inverse transformation.
As a final step, the restored image needs to be transformed into the original coordinate system by two times of coordinate inverse transformation.
The first one is polar coordinate inverse transformation. Pixel information is transformed from polar coordinates to the rectangular coordinate system of the plane π. The formula of the inverse transformation is shown in equation (25), and the schematic diagram of it is shown in figure 8, The second inverse transformation is to transform rectangular coordinates from the plane π to the imaging plane. The corresponding transformation matrix is the inverse matrix of the transformation matrix from the imaging plane to the plane π, which is shown in equation (21).
After completing all the above transformations, the final restoration result of the original star image with the elliptical  star streak can be obtained. Figure 9 shows the actual processing results after each step of the restoration algorithm.

Optimization of restoration results based on subdivid-
ing the polar coordinate. In the polar coordinate transformation, when the transformation is performed according to the conventional quantization unit (i.e. the rotation angle is quantized in the unit of 1 • and the rotation radius is quantized in the unit of 1 pixel), the restored star spot will suffer obvious tangential stretch deformation when transformed to the rectangular coordinate system, as shown in figure 8. The reason for this phenomenon is that the restored spot is not an ideal one. It is a point with a certain length and width. The width of the star spot depends on the quantization unit of the rotation angle. The larger the quantization unit is, the longer the width of the restored star spot is. The longer the width is, the longer the arc length of the star spot inverse-transformed to the rectangular coordinate system is. The larger the difference between the arc length and the radial width is, the more obvious the spot deformation is. Thereby, the optimal quantization unit of the rotation angle is selected to effectively solve the spot deformation, thus ensuring the restoration accuracy.
Firstly, the quantization unit of the rotation radius ∆r is set to 1 pixel. Then, the quantization unit of the rotation angle is ∆θ opt width of the recovered star spot in polar coordinate system is W r and length of it is L θ , where the length L θ corresponds to N rotation angle quantization units. In the rectangular coordinate system, the optimal quantization unit value of the tangential rotation angle ∆θ opt will be different if the star spot is located at a different radius r. Therefore, the radius is chosen to be half of the polar coordinate rotation radius, i.e. 1/4 length of the image sensor. The optimal rotation angle quantization unit corresponding to the above is used as the optimal rotation quantization unit for the transformation in this paper, so that the optimal rotation angle quantization unit value is obtained according to equation (26), In this paper, the radial width of the recovered star spot W r in polar coordinate system is three pixels, the length of the recovered star spot N is four quantization units. The simulated star sensor contains 2048 × 2048 pixels. Hence, r is taken as the middle value of 0-1024 pixels, i.e. 512 pixels. Therefore, in this paper, the optimal quantization unit value of the tangential rotation angle is 0.1 • . Figure 10 shows the processing results of polar coordinate transformation with the radial rotation radius quantized in 1 pixel and the rotation angle quantized in 0.1 • . As can be seen from figure 10, the transverse width of the restored star spot is significantly compressed compared with the conventional polar coordinate transformation method. And there is no tangential stretching deformation of the star spot after the polar coordinate inverse transformation, which improves the restoration accuracy.

Simulation analysis
The proposed algorithm is simulated and validated according to the typical parameters of the star sensor. The pixel array of the star sensor consists of 2048 × 2048 pixels. The focal length of the star sensor is 55 mm. The angle of FOV is 25 • . The exposure time of the star sensor is 0.1 s. The pixel gray value of the image at the time of exposure is determined according to equation (27). Typical star spot size is 3 × 3 pixels. The brightness level A of the star spot is set to 0.3 and the dispersion radius σ r is set to 0.35. N g is the added Gaussian noise with Based on the algorithm in this paper, simulation images at different angular velocities are tested and validated. Then, the results are compared with similar existing restoration algorithms. Figure 11 shows two representative results chosen from all experimental conditions. Figure 11(a) shows the result before and after restoration when the triaxial angular velocities are ω x = 5 rad s −1 , ω y = 5 rad s −1 , ω z = 5 rad s −1 , and the noise intensities are σ g = 0.00, σ g = 0.01, σ g = 0.015 respectively. Figure 11(b) shows result before and after restoration when the triaxial angular velocities are ω x = 5 rad s −1 , ω y = 5 rad s −1 , ω z = 5 rad s −1 , ω x = 5 rad s −1 , ω y = 5 rad s −1 , ω z = 10 rad s −1 , ω x = 5 rad s −1 , ω y = 5 rad s −1 , ω z = 15 rad s −1 respectively and the noise intensity is σ g = 0.01. From figure 11, it can be seen that the algorithm in this paper can effectively restore the star streaks at different compound angular velocities. Figure 12 shows the experimental results of the restoration accuracy with different triaxial angular velocities ⃗ ω ∈ Non-blind restoration Deep coupling of star tracker and MEMS-gyro [9] 0.16 pixel Multi-seed growth followed by RCRM restoration method [32] 0.1471 pixel Blind restoration Traditional motion blur restoration algorithm model [9] 0.36 pixel RL restoration method based on second-order vector extrapolation [23] 0.3225 pixel Real-time star tailing removal method based on fast blur kernel estimations [28] 0.2489 pixel Algorithm of this paper Noise intensity σg (0-0.03) 0.2370 pixel Angular velocity ⃗ ω (8-17 rad s −1 ) 0.2983 pixel Figure 13. Restoration accuracy at different noise intensities.
[8, 17] rad s −1 . It can be seen that the centroid localization accuracy of the restored images is about 0.19-0.39 pixels.
In the existing literature, there are only studies of linear motion-blurred restoration. Therefore, the accuracy of the proposed algorithm compared with similar existing methods, is shown in table 1. It can be seen that the accuracy is at the same level as existing star image restoration methods for straight star streaks.
Besides, the robustness of the proposed algorithm is tested and validated using images with different noise intensities. The noise is Gaussian noise with intensity σ g (0-0.03). The accuracy for each data group is shown in figure 13. It can be seen that the proposed algorithm is still able to restore effectively under different noise intensities.

Experimental validation
The algorithm of this paper is tested on actual images with laboratory star sensor and the corresponding test equipment. The experiment setup is shown in figure 14. The singlestar simulator is used to simulate a single-star target, and the triaxial turntable is used to simulate triaxial compound motion conditions. Figure 15(a) shows the actual image taken when the triaxial angular velocities are ω x = 1 rad s −1 , ω y = 1 rad s −1 , ω z = 15 rad s −1 respectively. Figure 15(b) shows the restoration image of the actual image with the proposed algorithm. It can be seen that the algorithm in this paper can restore the actual image effectively.
In the actual experimental test, due to the limited aperture of the single-star simulator (22 cm), if the angular velocity is too large, the star streak will quickly run out of the imaging plane. If the angular velocity is too small, the star streak is  approximately a straight star streak. Therefore, the experiment has some limitations. But it can still be able to validate the effectiveness of the proposed algorithm to some degree.

Conclusion
In this paper, a restoration algorithm for the star streak which is not a straight line is proposed. Through motion modeling analysis, it is found that the star streak is an elliptical arc under non-zero triaxial angular velocity. First, the elliptical star streak is transformed into a circular star streak by projective transformation. Then, the circular star streak is transformed into a straight star streak by polar coordinate transformation.
The original elliptical star streak is restored by restoration methods for straight star streaks and coordinate inverse transformation. At the same time, the algorithm is further optimized by subdividing the polar coordinates. The proposed algorithm is tested and validated by simulation and experimentation. The experimental results show that the algorithm is effective. The accuracy of the restored image is about 0.15-0.32 pixels with the noise intensity σ g being 0-0.03 and the triaxial angular velocity being ω x = 5 rad s −1 , ω y = 5 rad s −1 , ω z = 5 rad s −1 respectively. The accuracy of the restored image is about 0.19-0.39 pixel with the noise intensity σ g being 0.01 and the triaxial angular velocity being ⃗ ω ∈ [8, 17] rad s −1 respectively. The accuracy is at the same level as that of existing star image restoration methods for straight star streak.
This method can ensure the accuracy of star spot centroid extraction and altitude calculation under triaxial compound motion condition, and consequently improve the performance of star sensor.
In the future, more work will be done on this topic, especially in the direction of improving the restoring speed and increasing the robustness.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.