Star Centroiding Based on Fast Gaussian Fitting for Star Sensors

The most accurate star centroiding method for star sensors is the Gaussian fitting (GF) algorithm, because the intensity distribution of a star spot conforms to the Gaussian function, but the computational complexity of GF is too high for real-time applications. In this paper, we develop the fast Gaussian fitting method (FGF), which approximates the solution of the GF in a closed-form, thus significantly speeding up the GF algorithm. Based on the fast Gaussian fitting method, a novel star centroiding algorithm is proposed, which sequentially performs the FGF twice to calculate the star centroid: the first FGF step roughly calculates the Gaussian parameters of a star spot and the noise intensity of each pixel; subsequently the second FGF accurately calculates the star centroid utilizing the noise intensity provided in the first step. In this way, the proposed algorithm achieves both high accuracy and high efficiency. Both simulated star images and star sensor images are used to verify the performance of the algorithm. Experimental results show that the accuracy of the proposed algorithm is almost the same as the GF algorithm, higher than most existing centroiding algorithms, meanwhile, the proposed algorithm is about 15 times faster than the GF algorithm, making it suitable for real-time applications.


Introduction
Star sensors are indispensable components of spacecraft attitude measurement systems. They capture images of stars with an imaging chip, and determine the three-axis attitude according to the location of stars in the field of view (FOV) [1]. In recent years, with the rapid development of space exploration [2][3][4], higher attitude accuracy is being more and more required, making it imperative to improve the accuracy of star sensors. Many factors can affect the accuracy of a star sensor, including the size of FOV, sensitivity of the imaging chip, distortion of the lens, accuracy of star centroiding, etc. In this work, we focus on improving the star centroiding accuracy.
A typical method to get an accurate star centroid is to make the star sensor slightly defocused, making the size of a star spot not smaller than 3 × 3 pixels. The intensity of a de-focused star spot conforms to the "Airy disk", which can be approximated by the Gaussian function [5]. The task of a star centroiding algorithm is to estimate the centroid of the star spot.
Many algorithms have been developed for star centroiding, which can be divided into two classes: the gray scale centroid methods and the fitting methods. The gray scale centroid methods estimate the star centroid by calculating the center of gravity of the star spot, among which the most commonly used algorithm is the center of gravity algorithm (CG) [6]. The formula of the CG is: where i represents the i-th pixel; x c and y c represent the centroid of the star spot; I i represents the gray value of pixel i; x i and y i represent the coordinates of pixel i. The CG is simple and effective, but the accuracy of this algorithm is restricted by the S-curve error, which is introduced by discrete sampling of the star image [7]. A typical way to improve the accuracy of the CG is to insert a set of weight values in the formula of the CG, which is named the weighted center of gravity algorithm (WCG) [8]. The formula of the WCG is: where W i is the weight of the i-th pixel.
Theoretically, the weight values should be assigned by the Gaussian function [9], however, those values could not be easily obtained. In the squared gray centroid method [10], the weight W i is approximated by the intensity I i based on the fact that the intensity of a star spot conforms to the Gaussian function. The accuracy of this algorithm is significantly improved when the radius of the star spot is large enough, however, the accuracy decreases rapidly when the radius gets smaller.
The fitting methods estimate the star centroid by fitting the intensity of a star spot to the Gaussian function. As the star spot conforms to a 2-dimensional Gaussian function, the Gaussian fitting algorithm (GF) can in theory achieve the highest star centroiding accuracy [11][12][13]. However, since solving the Gaussian parameters is a nonlinear optimization problem, which is inevitably a multiple-step iteration process, the algorithm is quite sensitive to the initial parameters and is time consuming in practice. To address the shortcomings of the GF, the Gaussian analysis algorithm (GA) is proposed in [14], where it is assumed that there is no intensity noise in the star images and a closed-form solution of the Gaussian parameters is derived, significantly speeding up the GF, but the algorithm is sensitive to intensity noise. Another way to improve the GF is to divide the Gaussian fitting into multiple 1-dimensional Gaussian fittings [15]. These algorithms speed up the GF to some extent at the cost of decreasing the accuracy.
Besides the above two classes of star centroiding methods, some other new star centroiding methods were proposed recently. For instance, Flewelling in [16] presented a star centroiding method with information theoretic weighting. Although the method achieves quite high accuracy, the high computational complexity makes it unsuitable for real-time applications.
In this paper, we propose a star centroiding algorithm based on the fitting method. By transforming the objective function of the GF to a simpler form, the proposed algorithm can solve the Gaussian parameters without iteration, making it more efficient than the GF. On the other hand, the algorithm performs similarly as the GF in accuracy and noise robustness. This paper is organized as follows: in Section 2, we develop an algorithm to quickly solve the Gaussian parameters and propose a novel star centroiding algorithm. In Section 3, experiments are designed to verify the performance of the algorithm.

Imaging Characteristics of the Star Sensor
A star sensor is an imaging system that uses an optical system to image stars onto the array's photosensitive area. Its output is a series of numerical values that represent the scene intensity at a series of discrete locations [17]. The imaging model of star sensor can be formulated as: where I i , S i and N i represent the observed, source and noise intensity of pixel i, respectively. As described above, the source intensity of each pixel in the star spot can be expressed by a Gaussian function: where x i and y i represent the coordinates of pixel i; v represents the Gaussian parameters (A, x c , y c , σ x , σ y ); A represents the brightness level of the star; x c and y c represent the centroid of the star; σ x and σ y represent the standard deviation of the function.

Characteristics of Gaussian Fitting
As the most accurate star centroiding method, the GF estimates the star centroid by fitting the observed intensity of a star to the Gaussian function. The objective function of the GF is: where z i = I(x i , y i ) − S(x i , y i |v), and U is a set formed by the pixels in the star spot.
According to the theory of error, the star centroid estimated with Equation (5) is the optimum estimation in the maximum likelihood sense when the noise of each pixel conforms to a Gaussian distribution, which can be usually satisfied for star images. However, since Equation (5) is a nonlinear least squares problem, complex algorithms such as the Levenberg-Marquardt algorithm [18] or the trust-region-reflective algorithm [19] are needed to solve it. These algorithms usually start from an initial guess of the parameters and approach the optimal solution in a multiple iteration procedure [20], which is rather time-consuming. On the other hand, improper initial guess of the parameters may make the algorithms converge to local optimal parameters, leading to instability in the results [21].

Fast Gaussian Fitting
To overcome the shortcomings of the GF and preserve its accuracy, we develop the fast Gaussian fitting method in this section.
Modify the objective function of the GF by performing logarithm operation on both I and S in Equation (5) as follows: where g i = ln I(x i , y i ) − ln S(x i , y i |v). By substituting Equation (4) into Equation (6), it is easy to see that the optimization problem in Equation (6) is a linear least squares problem, which can be solved in a closed-form very efficiently. If Equation (6) shares the same solution with Equation (5) then it is a good acceleration method for the GF.
When there is no noise in images, i.e., N = 0, it is clear that Equation (6) is equivalent with Equation (5). However, in real star images, noise is always present and cannot be completely removed, i.e., N = 0. Therefore, the solution of Equation (6) is probably different from that of Equation (5).
By substituting Equation (3) into Equation (5) and Equation (6), the two objective functions can be converted to: where: This means Equation (8) can be seen as a weighted GF, and the weight of each pixel is determined by the intensity of source and noise of the pixel. For a certain noise intensity, a lower source intensity leads to a higher weight, which makes Equation (8) tend to be more sensitive to the pixels with lower source intensity. We illustrate this phenomenon with a simulated star spot in Figure 1. The intensity and weight of each pixel are shown in Figure 1a,b, respectively. Compared with the pixels in the center region, the pixels located at the edge of the star spot are much darker, but their weights are multiple times larger, which means the solution of Equation (8) is mainly determined by the edge pixels. Since the intensity of the edge pixels is easier to be influenced by noise, the center estimated with Equation (6) is more sensitive to noise than the GF, which is not what we expect. where: This means Equation (8) can be seen as a weighted GF, and the weight of each pixel is determined by the intensity of source and noise of the pixel. For a certain noise intensity, a lower source intensity leads to a higher weight, which makes Equation (8) tend to be more sensitive to the pixels with lower source intensity. We illustrate this phenomenon with a simulated star spot in Figure 1. The intensity and weight of each pixel are shown in Figure 1a,b, respectively. Compared with the pixels in the center region, the pixels located at the edge of the star spot are much darker, but their weights are multiple times larger, which means the solution of Equation (8) is mainly determined by the edge pixels. Since the intensity of the edge pixels is easier to be influenced by noise, the center estimated with Equation (6) is more sensitive to noise than the GF, which is not what we expect.    The intensity of the pixels in the center region is greater than that of the pixels in the edge region in (a). The weights of pixels in the center region are much less than that of pixels in the edge region in (b). All the values in (c) are very close to 1.  The intensity of the pixels in the center region is greater than that of the pixels in the edge region in (a). The weights of pixels in the center region are much less than that of pixels in the edge region in (b). All the values in (c) are very close to 1.
To address this problem, we need to eliminate the nonuniformity of the weights. g i in Equation (6) can be transformed into: , and define |r i | as the signal-to-noise ratio (SNR) of the i-th pixel. ϕ S i N i can be transformed into: When r i of pixel i ranges from −1 to 1, i.e., SNR of pixel i is less than 1, the pixel is submerged in the noise and is usually excluded from the star spot area, thus there is no need to analyze this case. When SNR of pixel i is greater than 1, the value of ϕ(r) converges to 1 with the increase of the SNR, as shown in Figure 2.
When i r of pixel i ranges from 1 − to 1, i.e., SNR of pixel i is less than 1, the pixel is submerged in the noise and is usually excluded from the star spot area, thus there is no need to analyze this case. When SNR of pixel i is greater than 1, the value of ( ) r ϕ converges to 1 with the increase of the SNR, as shown in Figure 2. Considering the values of ( ) r ϕ of pixels in the star spot are close to 1, we modify Equation (6) by adding weight i I in the equation as follows: where U is a set formed by the pixels whose ( ) r ϕ are close to 1.
Equation (12) can also be seen as a weighted GF with the weight of each pixel determined by Figure 1c illustrates the values of ( ) r ϕ for the star spot in Figure 1a. Different from the weights in Figure 1b, all the values of ( ) r ϕ in Figure 1c are very close to 1, indicating that the solution of Equation (12) can be approximately equal to Equation (5), and is more robust to noise than that of Equation (6). Considering the values of ϕ(r) of pixels in the star spot are close to 1, we modify Equation (6) by adding weight I i in the equation as follows: where U is a set formed by the pixels whose ϕ(r) are close to 1. Equation (12) can also be seen as a weighted GF with the weight of each pixel determined by ϕ(r i ). Figure 1c illustrates the values of ϕ(r) for the star spot in Figure 1a. Different from the weights in Figure 1b, all the values of ϕ(r) in Figure 1c are very close to 1, indicating that the solution of Equation (12) can be approximately equal to Equation (5), and is more robust to noise than that of Equation (6). It can be seen from Figure 2, when the SNR of a pixel is greater than a certain value, the value of ϕ(r) can be limited in any given neighborhood of 1. For example, when SNR > 3, ϕ(r) ∈ (0.8109, 1.1507); when SNR > 10, ϕ(r) ∈ (0.9482, 1.0484). Therefore, we can select the pixels with SNR greater than a certain value to ensure that ϕ(r) of these pixels is close enough to 1.
Moreover, Equation (12) is also a linear least squares problem, which can be solved efficiently in closed-form, making it a good acceleration method for the GF. We name the method as the fast Gaussian fitting method. The process of finding the optimal solution of Equation (12) is described below. By substituting Equation (4) to Equation (12), Equation (12) is converted into: For the convenience of calculations, m, n, p, q and k are defined follows: By substituting Equation (14) to Equation (13), the objective function is converted into: For the parameters m, n, p, q and k, Equation (15) is a least squares problem. The solution can be found by solving the following equation: By using the derivative operation, we can find the critical point of each parameter, and solve the optimal parameters m, n, p, q and k. The corresponding parameters (A, x c , y c , σ x , σ y ) can be solved as follows:

Star Centroiding Based on Fast Gaussian Fitting
As described above, the objective function of the fast Gaussian fitting method can be approximately equal to that of the GF under the precondition: ϕ(r) of each pixel involved in the calculation is close enough to 1. In star images, due to the uncertainty of the noise and star intensity, not all pixels satisfy the precondition, and it does not provide how to select the pixels with SNR greater than a certain value. Therefore, the fast Gaussian fitting method cannot be used as a complete star centroiding method. In this section, we design a method to select the pixels that satisfy the precondition by thresholding the SNR of pixels, and propose a star centroiding algorithm based on the fast Gaussian fitting method. The proposed star centroiding algorithm is summarized as follows: when a star spot is extracted, two steps are designed to solve the star spot centroid. In the first step, we get a rough solution of the Gaussian parameters of the star spot with the FGF, then use the solution to estimate the noise and SNR of each pixel in the star spot. In this way, the SNR of each pixel in the star spot is obtained, and the pixels that satisfy the precondition in the FGF can be found by thresholding the SNR. In the second step, we select all the pixels that satisfy the precondition and use the FGF to calculate the star centroid. Compared with the FGF in the first step, more pixels are involved in the second step, and the solution is more accurate. The working procedure is shown in Figure 3.
approximately equal to that of the GF under the precondition: ( ) r ϕ of each pixel involved in the calculation is close enough to 1. In star images, due to the uncertainty of the noise and star intensity, not all pixels satisfy the precondition, and it does not provide how to select the pixels with SNR greater than a certain value. Therefore, the fast Gaussian fitting method cannot be used as a complete star centroiding method. In this section, we design a method to select the pixels that satisfy the precondition by thresholding the SNR of pixels, and propose a star centroiding algorithm based on the fast Gaussian fitting method.
The proposed star centroiding algorithm is summarized as follows: when a star spot is extracted, two steps are designed to solve the star spot centroid. In the first step, we get a rough solution of the Gaussian parameters of the star spot with the FGF, then use the solution to estimate the noise and SNR of each pixel in the star spot. In this way, the SNR of each pixel in the star spot is obtained, and the pixels that satisfy the precondition in the FGF can be found by thresholding the SNR. In the second step, we select all the pixels that satisfy the precondition and use the FGF to calculate the star centroid. Compared with the FGF in the first step, more pixels are involved in the second step, and the solution is more accurate. The working procedure is shown in Figure 3.  (1) Star spot extraction: When a star image is obtained, the connected components are extracted to locate the star spots in the image. The number of pixels in each star spot should be no less than 5, otherwise, the star spot will be considered as a false star and removed. In addition, the pixels with full saturated gray value in the star spot are excluded. (2) The first step: The task of this step is to estimate the SNR of each pixel. The detailed process is as follows: (1) Star spot extraction: When a star image is obtained, the connected components are extracted to locate the star spots in the image. The number of pixels in each star spot should be no less than 5, otherwise, the star spot will be considered as a false star and removed. In addition, the pixels with full saturated gray value in the star spot are excluded. (2) The first step: The task of this step is to estimate the SNR of each pixel. The detailed process is as follows: (i) Initial selection of pixels: Although the noise of each pixel is unknown, the pixels with larger gray values usually have higher SNR, thus their corresponding ϕ(r) are closer to 1. According to this idea, the five pixels with the maximum gray values in the star spot are selected as the initial pixels of the set U.

(ii)
Fast Gaussian fitting: The fast Gaussian fitting method described above is adopted to find the solution of the Gaussian parameters (A , x c , y c , σ x , σ y ). With these parameters, we can establish the star intensity function, which is: (iii) Pixel-wise SNR estimation: According to Equation (3) and Equation (18), the noise intensity of each pixel in the star spot can be calculated by Equation (19), and the SNR of each pixel can be estimated by Equation (20): (3) The second step: Since the SNR of each pixel has been estimated in the previous step, we can solve the star centroid as follows: (i) Reselect pixels based on SNR: In this step, we will update the pixels in the set U by thresholding the SNR. The pixels in the star spot are filtered by judging whether the SNR of the pixel is greater than a threshold and those pixels with SNR greater than the threshold are collected to form the set U. T is used as the value of the threshold: (ii) Fast Gaussian fitting: The fast Gaussian fitting method is adopted to find the optimal solution of the Gaussian parameters again. Comparing with the first step in the noise estimation step, more pixels are involved in this step, and the results will be more accurate.

Results and Discussion
To evaluate the performance of the proposed algorithm, we test the algorithm on a set of simulated star images with different characteristics, and compare with the existing star centroiding methods. In the following sub-sections, we describe generation of the simulated star images, selection of parameters, and experiments in detail.

Star Images Generation
In this work, the star images are generated to simulate the images obtained by a star sensor. Each star image is of 31 × 31 pixels, and includes only one star spot.
In the simulated star images, each pixel is generated according to Equation (4), as follows: where x i and y i represent the i-th pixel in the image; x c and y c represent the centroid of the star spot in pixels; σ x = σ y = σ r represents the standard deviation of the generated function of the star spot, and is defined as the Gaussian radius here; A represents the brightness level here; N i represents the noise intensity of pixel i. The value of N i is a random value which obeys Gaussian-distributed with 0 mean and σ g standard deviation.

Parameter Selection
The threshold T in the fast Gaussian fitting method is a key parameter of the proposed algorithm. We select the value of the parameter as follows: The threshold T is used to select the pixels involved in the star centroiding step. Two aspects should be taken into account to select the value of T: (i) The number of pixels in set U must be no less than 5; (ii) ϕ(r) of each pixel selected should be close enough to 1.
If the value of T is very large, ϕ(r) of each pixel in set U is closer to 1, but the number of pixels in U is reduced. On the other hand, if the value of T is very small, the number of the pixels in U increases, but ϕ(r) of some pixels in U may be far away from 1, which can reduce the accuracy of the fast Gaussian fitting method. Considering the diversity of the star images, we set T = 3.

Accuracy Experiments
The proposed algorithm is verified by a series of star images with different centroid, Gaussian radius, noise level and brightness level. The accuracy of the algorithm is analyzed by the star centroiding error, which is defined as the Euclidean distance between the calculated centroid (x c ,ŷ c ) and the true centroid (x c , y c ): Several existing star centroiding methods are adopted to compare with the proposed algorithm. The center of gravity algorithm (CG) is considered to be the least time-consuming method, and the squared gray centroid method is the most common weighted center of gravity algorithm (WCG), which is an improvement of the CG. The Gaussian fitting algorithm (GF) is the most accuracy method in theory, and the Gaussian analytic method (GA) is the latest improved Gaussian Fitting algorithm, which is claimed to be the most accurate and efficient algorithm when the size of the star spot radius is not bigger than 5 × 5. All of the experiments are carried on the MATLAB 2016b software platform running on an Intel Core i5-7500T 2.7 GHz processor.

Experiment with the Star Spot at Different Locations
In this experiment, the accuracy of the algorithm is verified by moving the centroid of the star spot from one side of a pixel to the other. Considering the symmetry of the star spot, the star spot is set to move along the x-axis at the center of the pixel, which is located at the center of the image, i.e., x o ranges from 15 to 16, and y o = 15.5. The Gaussian radius of the star spot is 0.5, the brightness level is 0.7, and the standard deviation of the Gaussian noise is 0.001, which is equivalent to σ = 4.096 for a 12-bit AD CMOS. For each location of the star spot, 10,000 images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 4. This phenomenon is caused by the S-curve error, which is an artefact of the sampling theorem due to undersampling [7,22]. In this experiment, the S-curve error is clearly displayed and is found only in the gray scale centroid methods. There is no significant S-curve error in the centroiding error of the fitting methods. There is a sine curve in the pixel for the CG and WCG, and the error of the WCG is greater than that of the CG.

Experiment with Different Gaussian Radius
In this experiment, the accuracy of the algorithm is tested by changing the Gaussian radius of the star spot in the simulated images. The brightness level of each star spot is fixed at 0.7 and the Gaussian radius varies from 0.5 to 2. Gaussian noise with 0.001 g σ = is also added to each image.
The centroid of the star spot is uniformly distributed in central pixel of the star image. For each radius, 10,000 simulated star images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 5.
centroiding error There are three special points (15, 15.5), (15.5, 15.5) and (16, 15.5) in this experiment. When the star spot is located in any of the three points, the centroiding error of the CG gets its local minimum value. When the star spot moves along the x-axis in the pixel, the centroiding error curve of the CG is a sine-like curve. There is a similar phenomenon for the WCG, whose error is greater than that of the CG. There is no obvious such phenomenon for the GF, the GA and the proposed star centroiding algorithm, whose error curves are stable when the star spot moves from one side to another. This phenomenon is caused by the S-curve error, which is an artefact of the sampling theorem due to undersampling [7,22]. In this experiment, the S-curve error is clearly displayed and is found only in the gray scale centroid methods. There is no significant S-curve error in the centroiding error of the fitting methods.

Experiment with Different Gaussian Radius
In this experiment, the accuracy of the algorithm is tested by changing the Gaussian radius of the star spot in the simulated images. The brightness level of each star spot is fixed at 0.7 and the Gaussian radius varies from 0.5 to 2. Gaussian noise with σ g = 0.001 is also added to each image. The centroid of the star spot is uniformly distributed in central pixel of the star image. For each radius, 10,000 simulated star images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 5. There is a sine curve in the pixel for the CG and WCG, and the error of the WCG is greater than that of the CG.

Experiment with Different Gaussian Radius
In this experiment, the accuracy of the algorithm is tested by changing the Gaussian radius of the star spot in the simulated images. The brightness level of each star spot is fixed at 0.7 and the Gaussian radius varies from 0.5 to 2. Gaussian noise with 0.001 g σ = is also added to each image.
The centroid of the star spot is uniformly distributed in central pixel of the star image. For each radius, 10,000 simulated star images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 5. Due to the existence of the S-curve error, the centroiding errors of the CG and the WCG are much greater than that of other algorithms when the Gaussian radius is 0.5. As the S-curve error decreases with the increase of the Gaussian radius [7], the centroiding errors of the CG and the WCG first decrease and then tend to be stable. When the error curve is stable, the random noise is the main component of the centroiding error, and the centroiding error of the WCG is less than that of the CG.
As the centroiding errors of the GA, the GF and the proposed algorithm are mainly composed of the random error, the error curves of those algorithms are more stable as shown in Figure 5. When the Gaussian radius is less than 0.6, the random error increases with the decrease of the pixels in the star spot, which leads to the increase of the centroiding error of the GA, the GF and the proposed algorithm.
The GA performs better than the CG and the WCG when the radius of the star spot is relatively small, but its centroiding error is greater than other algorithms when the Gaussian radius of the star spot is greater than 1.4. The centroiding error of the proposed algorithm is almost the same as that of the GF under all the radius of the star spot, less than that of other algorithms.

Experiment with Different Noise Level
In this experiment, the accuracy of the algorithm is verified by changing the level of the noise added in the star images. The standard deviation of the noise ranges from 0 to 0.03. The Gaussian radius of the star spot is uniformly distributed between 0.5 and 1.2. The brightness level is 0.7. The centroid of the star spot is uniformly distributed within the central pixel of the star image. For each level of the noise, 10,000 images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 6. star spot, which leads to the increase of the centroiding error of the GA, the GF and the proposed algorithm.
The GA performs better than the CG and the WCG when the radius of the star spot is relatively small, but its centroiding error is greater than other algorithms when the Gaussian radius of the star spot is greater than 1.4. The centroiding error of the proposed algorithm is almost the same as that of the GF under all the radius of the star spot, less than that of other algorithms.

Experiment with Different Noise Level
In this experiment, the accuracy of the algorithm is verified by changing the level of the noise added in the star images. The standard deviation of the noise ranges from 0 to 0.03. The Gaussian radius of the star spot is uniformly distributed between 0.5 and 1.2. The brightness level is 0.7. The centroid of the star spot is uniformly distributed within the central pixel of the star image. For each level of the noise, 10,000 images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 6. As seen in Figure 6, with the increase of the noise level, the centroiding error of each algorithm increases gradually. When the noise level g σ is greater than 0.006, the centroiding errors of the CG and the GA are similar, both greater than that of other algorithms, indicating that the two algorithms are less robust to noise. The WCG performs better than the CG and the GA when the noise level g σ is relatively high, but its centroiding error is the worst when the noise level g σ is less than 0.004.
The centroiding error of the proposed algorithm is almost the same as that of the GF, less than other algorithms under all levels of the noise. It indicates that the robustness to noise of the proposed algorithm is almost the same as that of the GF, better than other star centroiding algorithms. As seen in Figure 6, with the increase of the noise level, the centroiding error of each algorithm increases gradually. When the noise level σ g is greater than 0.006, the centroiding errors of the CG and the GA are similar, both greater than that of other algorithms, indicating that the two algorithms are less robust to noise. The WCG performs better than the CG and the GA when the noise level σ g is relatively high, but its centroiding error is the worst when the noise level σ g is less than 0.004.
The centroiding error of the proposed algorithm is almost the same as that of the GF, less than other algorithms under all levels of the noise. It indicates that the robustness to noise of the proposed algorithm is almost the same as that of the GF, better than other star centroiding algorithms.

Experiment with Different Brightness Level
In this experiment, the accuracy of the algorithm is verified by changing the brightness level of the stars in the images. The brightness level ranges from 0.15 to 2 (The saturation brightness of each pixel is 1). The Gaussian radius of the star spot is 0.7. The centroid of the star spot is uniformly distributed within the central pixel of the star image. The Gaussian noise with σ g = 0.001 is added to each image. For each level of the brightness, 10,000 images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 7.
the stars in the images. The brightness level ranges from 0.15 to 2 (The saturation brightness of each pixel is 1). The Gaussian radius of the star spot is 0.7. The centroid of the star spot is uniformly distributed within the central pixel of the star image. The Gaussian noise with 0.001 g σ = is added to each image. For each level of the brightness, 10,000 images are generated, and the mean centroiding error of the 10,000 images is used as the centroiding error. The error curves are shown in Figure 7. With the increase of the brightness, there are two main factors affecting the centroiding accuracy: (i) When a pixel in the star spot is unsaturated, the SNR of the pixel increases with the star getting brighter, helping to improve the accuracy of the algorithm. (ii) When a pixel is saturated, the gray value of the pixel is truncated, which makes it different from the true value, leading to the increase of the centroiding error.
In Figure 7, factor (i) is the main factor when 1 A ≤ , the centroiding error of all the algorithms in this phase reduce because of the increase of SNR. When A continues to get larger, factor (ii) becomes the main factor, the centroiding errors of the CG, the WCG and the GA increase seriously because of the truncation of the saturated pixels. In contrast, the centroiding error of the proposed algorithm and the GF keep stable in this phase. The robustness of the two algorithms results from that the pixels with saturated gray value in the star spot are excluded, as described in Section 2.4.

Efficiency Experiment
The efficiency of the algorithm is also tested and compared with the other algorithms. The total time consumption of each algorithm on 10,000 star images is counted and compared, and listed in Table 1. The brightness level is 0.7 and the Gaussian radius of each star spot is uniformly distributed between 0.5 and 1.2, and the position of the star spot is uniformly distributed within the central pixel of the star image and the standard deviation of the noise is randomly obtained from 0 to 0.02.  With the increase of the brightness, there are two main factors affecting the centroiding accuracy: (i) When a pixel in the star spot is unsaturated, the SNR of the pixel increases with the star getting brighter, helping to improve the accuracy of the algorithm. (ii) When a pixel is saturated, the gray value of the pixel is truncated, which makes it different from the true value, leading to the increase of the centroiding error.
In Figure 7, factor (i) is the main factor when A ≤ 1, the centroiding error of all the algorithms in this phase reduce because of the increase of SNR. When A continues to get larger, factor (ii) becomes the main factor, the centroiding errors of the CG, the WCG and the GA increase seriously because of the truncation of the saturated pixels. In contrast, the centroiding error of the proposed algorithm and the GF keep stable in this phase. The robustness of the two algorithms results from that the pixels with saturated gray value in the star spot are excluded, as described in Section 2.4.

Efficiency Experiment
The efficiency of the algorithm is also tested and compared with the other algorithms. The total time consumption of each algorithm on 10,000 star images is counted and compared, and listed in Table 1. The brightness level is 0.7 and the Gaussian radius of each star spot is uniformly distributed between 0.5 and 1.2, and the position of the star spot is uniformly distributed within the central pixel of the star image and the standard deviation of the noise is randomly obtained from 0 to 0.02. There is no doubt that the CG is the fastest. The GF, which takes about 44 times as much time as the CG is the slowest. The WCG and the GA are slightly slower than the CG. The proposed algorithm takes about 3 times the amount of time of the CG, and is about 15 times faster than the GF.

Experiment with Star Sensor Imagery
To further verify the performance of the proposed algorithm, both the CG and the proposed algorithm are implemented and tested on an existing star sensor. Because the true value of the star centroid cannot be obtained, the accuracy of the algorithm is difficult to evaluate. In this section, repeatability is used to evaluate the performance of those two algorithms, and an experiment is performed utilizing a starlight simulator in laboratory.
The star sensor indoor testing system contains a starlight simulator, a two-axis turntable and a tested star sensor. The starlight simulator provides highly collimated starlight with different magnitude, and the turntable provides highly precision test position by rotating its internal frame and external frame. The tested star sensor is mounted on the internal frame of the turntable with its boresight pointing at the starlight simulator. The experiment is carried in a dark room.
By rotating the internal frame and the external frame of the turntable, the star is imaged at different positions of the image plane. For each position, 100 star images are captured and processed by both the two algorithms.
As the real position of the star in the image plane cannot be obtained, the performance of the two algorithms is compared according to the standard deviation of the centroids at each position. The algorithm with smaller standard deviation has better repeatability. The comparison results of all the positions are shown in Figure 8. There is no doubt that the CG is the fastest. The GF, which takes about 44 times as much time as the CG is the slowest. The WCG and the GA are slightly slower than the CG. The proposed algorithm takes about 3 times the amount of time of the CG, and is about 15 times faster than the GF.

Experiment with Star Sensor Imagery
To further verify the performance of the proposed algorithm, both the CG and the proposed algorithm are implemented and tested on an existing star sensor. Because the true value of the star centroid cannot be obtained, the accuracy of the algorithm is difficult to evaluate. In this section, repeatability is used to evaluate the performance of those two algorithms, and an experiment is performed utilizing a starlight simulator in laboratory.
The star sensor indoor testing system contains a starlight simulator, a two-axis turntable and a tested star sensor. The starlight simulator provides highly collimated starlight with different magnitude, and the turntable provides highly precision test position by rotating its internal frame and external frame. The tested star sensor is mounted on the internal frame of the turntable with its boresight pointing at the starlight simulator. The experiment is carried in a dark room.
By rotating the internal frame and the external frame of the turntable, the star is imaged at different positions of the image plane. For each position, 100 star images are captured and processed by both the two algorithms.
As the real position of the star in the image plane cannot be obtained, the performance of the two algorithms is compared according to the standard deviation of the centroids at each position. The algorithm with smaller standard deviation has better repeatability. The comparison results of all the positions are shown in Figure 8. The points in Figure 8a represent the positions of the stars captured in experiment. The points marked with red asterisks indicate that the proposed algorithm performs better than the CG there, The points in Figure 8a represent the positions of the stars captured in experiment. The points marked with red asterisks indicate that the proposed algorithm performs better than the CG there, while the points marked with blue dots indicate that the CG performs better. It can be seen that the proposed algorithm performs better at most of the positions. The positions that the CG performs better locate mainly near the edge of the FOV. One typical sTabletar spot at these positions is shown in Figure 8c. Due to the distortion of the optical system, the intensity of the star spot does not conform to Gaussian function strictly, making the proposed algorithm perform slightly worse than CG sometimes. Figure 8b is a star spot located near the center of the FOV, the intensity of which conforms to the Gaussian function better. Table 2 shows the standard deviation of the star spots in Figure 8b,c.

Conclusions
In order to improve the accuracy of a star sensor, this work focuses on improving the accuracy of the sensor. The Gaussian fitting algorithm is a star centroiding method with high accuracy and low efficiency, and the star centroiding algorithm we proposed aims to improve the efficiency of the Gaussian fitting algorithm without reducing the accuracy. For star images with different centroid, radius, brightness level and noise level, the proposed algorithm performs better than the existing star centroiding methods in terms of accuracy, robustness and efficiency. The proposed algorithm is also evaluated on an existing star sensor, and the results indicate that its repeatability is better than the center of gravity algorithm in most cases. The accuracy and efficiency of this algorithm make it suitable for real-time applications.