POSITION ESTIMATION USING IMAGE DERIVATIVE

This paper describes an image processing algorithm to process Moon and/or Earth images. The theory presented is based on the fact that Moon hard edge points are characterized by the highest values of the image derivative. Outliers are eliminated by two sequential ﬁlters. Moon center and radius are then estimated by nonlinear least-squares using circular sigmoid functions. The proposed image processing has been applied and validated using real and synthetic Moon images.


INTRODUCTION
This paper completes a study, previously performed [1], to process cropped Moon and/or Earth images. The theory presented, which has been applied to real and synthetic Moon images, is based on the fact that the edge points are characterized by the highest values of the image derivative. The content of this paper follows the general flowchart given in Fig. 1 which summarizes the main processes. These steps are explained in details in the following sections. The proposed algorithm has the following steps: • The geometrical illumination conditions and the expected size of the observed Moon are first computed. These data allow to understand if the Moon is new, crescent, gibbous, or full.
• The image differentiation is obtained with a fourth order Richardson extrapolation. An analytical surface example is provided quantifying the accuracy gain obtained using Richardson extrapolation.
• The pixels of the gradient image thus obtained, are sorted in descending mode. This means that the first set of pixels have high probability of being edge points. Outliers are then eliminated by two sequential filters. The first filter validates an edge pixel based on eigenanalysis of a small gradient box around the pixel and by investigation of the same box in the original image. The second filter is an improved version of the RANSAC algorithm.
• The edge point selected are then used to obtain a first estimate of the Moon center and radius. Taubin best fitting [6] (for circles) and the Fitzgibbons algorithm (for ellipses) are used. For the current application, the eccentricity of the Earth is so small that the ellipse angular orientation is a parameter ill defined to estimate (using noisy edge points). In addition, the presence of the atmosphere and high-altitude clouds make the problem even more inaccurate for Earth images. Therefore the orientation of the Earth ellipse is assumed known from the attitude of the camera.
• The final estimation of Moon center and radius is performed by non-linear least-squares using a circular sigmoid function to describe the Moon edge transition.

ILLUMINATION AND IMAGE PROCESSING PARAMETERS
The image processing software require the computation of some internal parameters. With reference to Fig. 2, these parameters depend on the Moon size on the imager, p r (cannot be smaller than a given value) and on the illumination geometry, p i (to define the separation between new and crescent Moon and between gibbous and full Moon). In this section, based on the geometry (Moon and Sun positions, J2000) and on a rough estimate of Spacecraft position (J2000), the conditions when the Moon is geometrically considered new, full, gibbous, or crescent are here derived. The observed Moon radius (in pixel) in the imager, p r , as well as the illumination size along the axis of symmetry, p i , dictate the values of important internal driving parameters used in the image processing code.
To estimate the full or new Moon cases, we use the following input data. Letr s be Sun-to-Moon rays direction (unit-vector), r o be (estimated) Spacecraft-to-Moon vector, and R m the Moon radius. From these data we can compute the illumination angle is then computed as from which ϑ is derived. The angle α m can be derived from the equation Figure 3. Full/New Moon geometry Full or new Moon while the value of the critical distance, r cr , can be obtained by Therefore, if |r o | < r cr and ϑ < π/2 then p i = 2p r (Full Moon) ϑ > π/2 then p i = 0 (New Moon) (5)

Gibbous Moon
Gibbous Moon will occur when |r o | > r cr and ϑ < π 2 .
When these conditions are satisfied then (see Fig. 4) the p i illumination parameter can be computed using the following equations and and, therefore, Crescent Moon will occur when |r o | > r cr and ϑ > π 2 . (11) When these conditions are satisfied then (see Fig. 5) the p i illumination parameter can be computed using the following equations and and, therefore,

GRADIENT-BASED IMAGE PROCESSING ALGORITHM
The gradient-based image processing approach is based on the fact that the part of the image with strongest variation of gray-tones is the Moon observable hard edge. Pixels with highest derivatives are the most likely to belong to the hard edge. Some of the selected bright gradient pixels do not belong to the Moon hard edge (outliers). These outliers are eliminated prior to using a best fitting algorithm to determine the sought geometrical properties.

Image differentiation
The gradient of a discrete function (the imager) is computed by numerical differentiation. In our case, the gradient is computed using four point central difference with a single Richardson extrapolation.
The row and column four point central partial differences computed at pixel [r, c] are These derivatives are accurate with order h 4 , where h is the pixel dimension. The image gradient is then computed as g (r, c) = g 2 r (r, c) + g 2 c (r, c) A single Richardson extrapolation requires computing Eq. (17) one time with step 2 h (getting g 2h ) and one time with step h (getting g h ). Therefore, the true derivative can be written as Equating these two equation we obtain and therefore An example of application of a single Richardson extrapolation is done using the six-hump Camel back function   The image gradient pixels are sorted in descending mode. The pixels with highest gradient values (edge pixels) are, therefore, concentrated as the first sorted elements. Unfortunately, some high gradient pixels (outliers) also appear on the Moon body (because of the illumination and local topology) or external (e.g., some bright stars or visible planets). Therefore, the outliers must be eliminated. Two sequential approaches have been implemented in order to eliminate the outliers.

Box-based outliers Identification
The first approach to remove outliers is based on local analysis around the pixels being considered as potential edge pixels (pixels with the highest gradient values). With reference to Fig. 8, this approach is based on using a mask size (the figure shows a 7 × 7 mask) whose size depend on the expected Moon dimension in pixels, p r .
In Fig. 8 the central pixel (in blue) is the one being considered as potential pixel on the Moon edge. Using all the pixels of the gradient mask as lumped masses, the eigenvalues and eigenvectors of the inertia matrix are computed. The eigenvalues ratio (min/max) must be lower than a given threshold (e.g., 0.3), the eigenvector associated to the minimum eigenvalue (red line) in inclined by an angle γ with respect to the row axis direction. In the figure, the pixels marked as green should both have high derivative values. The pixels marked as red, on the other hand should have very different gray tones.

RANSAC for circle and ellipse
A standard algorithm to eliminate outliers is called RANSAC which is an abbreviation for "RANdom SAmple Consensus". The RANSAC algorithm [2] is an algorithm for robust fitting of models in the presence of many data outliers. The algorithm is very simple and can be applied ti lines, circle, ellipse, etc. Given a fitting problem with n data points and unknown m parameters x, estimate the parameters. For example, for lines m = 2, circles m = 3, and for ellipse m = 5. Starting with N min = n, the classic RANSAC algorithm perform N test times the following steps: 1. selects m data points at random; 2. estimates the m parameter, meaning to find the curve passing through these m points; 3. counts how many data items are located at a distance greater than a minimum value d min from the curve. They are N k .
4. if N k < N min then set N min = N k . One of the limitations of the RANSAC algorithm is the random procedure in which the points are selected. The random selection approach may persistently select outliers. To improve RANSAC, a technique (originally developed for Pyramid Star-Identification [3]) to select always different triads with indices continuously changing, has been adopted.

Best fitting for Circle and Ellipse
Once the pixels in the image have been selected, the algorithm must determine the geometric characteristics of the observed body, that is, radius and center for the Moon, considered a sphere, and axes and center for the Earth, considered an ellipsoid of revolution. Feature recognition for conic curves is based essentially on best fitting, which are typically grouped in two categories: Geometric fit tries to minimize the geometric distance between the data points and the conic. In general this method is considered the most accurate, however it does not admit closed solution and requires the use of iterative algorithms.
Algebraic fit tries to minimize the "algebraic distance", i.e. the implicit conic equation F (x, y) = 0. Many methods based on this approach have been developed and they differ in accuracy, stability and complexity. They generally give results less accurate than the geometric fit, but on the other hand have closed form solutions and therefore can be implemented in a much simpler way.

Taubin for circles
Best fitting algorithms based on geometric fit are iterative methods requiring an initial guess, and their reliability and accuracy is ultimately solely dependent in the iterative scheme implemented. Typically, the Levenberg-Marquardt scheme is considered the most accurate. However, given the need to keep the computational load on-board at a minimum, algorithms based on algebraic fit seem a better choice, even if their accuracy is slightly worse. Among these, the Taubin [6] method is used and described in further detail.
Given the circle implicit equation as which is equivalent to minimize subject to the constraint where the bar sign indicates simple average of the variables.
One important property of this method is that the results are independent of the reference frame. By reformulating this problem in matrix form, it becomes apparent how it is equivalent to a generalized eigenvalue problem. Indeed, define then the best fitting problem can be rewritten as where Therefore, the optimization problem is equivalent to the following generalized eigenvalue problem: Since for each solution, then to minimize A T M A the smallest positive eigenvalue should be chosen. Because the matrix N has rank three, such eigenvalue can be found by explicitly solving the cubic characteristic equation without the necessity to implement SVD or other more sophisticated algorithms.

Fitzgibbons for ellipse
The problem of best fitting for an ellipse is sensibly more complicated and somewhat less studied than the circle. Firstly, compared to the three parameters necessary to describe the circle (x, y, R), a generic ellipse is defined by 5 parameters (x, y, a, b, θ). Moreover, especially for very small values of the eccentricity, the inclination is very hard to discern and the error associated with it is several orders of magnitude greater than for the other parameters.
Because the Earth itself, when considered an ellipsoid of revolution, has very small eccentricity, its planar projections are also almost circular, with ratio between the axes varying between 1 and 0.9966 depending on the latitude of the observer. We thus expect the highest errors in inclination estimation.
While geometric fit methods exist for the ellipse, they tend to not be very accurate and efficient due to the complicated formulation required. An non-linear Least Square approach, using the Jacobian of an ellipse in canonical form has been attempted; however, this method has proven itself not accurate enough for the requirements, and extremely slow to converge, even when advanced algorithms, like Levenberg-Marquardt are used. On the other hand, several studies ( [4], [5]) have explored and improved the algebraic fit approach, with promising results. In particular, we have selected the Improved Fitzgibbon method, which reduces to solving a 3 × 3 eigenvalue problem.
For this equation to represent an ellipse, it must be b 2 − 4ac < 0, but with proper rescaling of the coefficient, which does not change the final result, it can be rewritten as b 2 − 4ac = 1. Therefore the problem can be rewritten in matrix form as:  Solving this optimization problem leads to Sa = λCa subject to a T Ca = 1, and the solution is the eigenvector corresponding to the minimum positive eigenvalue. However, solving the problem in this form leads to numerical instability. Instead, it is possible to decouple the quadratic part of D from the linear part, to obtain a simplified formulation that, together with being stable, has the added advantage of being of dimension 3 and thus the subsequent eigenanalysis can be performed without using SVD or similar. Therefore, by defining and where The final system of equation assumes the form: and This version of the Fitzgibbon algorithm has been implemented with an explicit cubic solver for the characteristic equation.

NUMERICAL TESTS
To support the capability of the developed image processing code, three numerical examples of partially observed Moon are here presented. All these tests have been performed with no information about time, camera data, and observer, Moon, and Sun positions. The image derivative is obtained using the 4-point central approach with no Richardson extrapolation. These information clearly help to define the best values of the internal parameters. Since this help is missing, the internal parameters have been adopted with same values in all three tests. These values are: Example #1: Real image, four-times cropped The original image of example #1, shown in Fig. 10, was intentionally cropped from a real image (taken from Earth) to validate the capability of the software to face a multiple cropped Moon images. Four times cropped is the also the maximum number of cropped segments possible. The four figures of Fig. 10 show the four main phases of the image processing. Once the image differentiation has been computed, the first 1,490 points with highest derivatives have been selected. This number differs from N = 1500 because the image is small and the number of points selected is also a function of the total number of pixels (N = 1500 is also the maximum number). The box-based outliers elimination approach discarded 375 points based on the results shown in Fig. 12 while RANSAC eliminated additional 5 points. The resulting points are shown in the bottom right figure of Fig. 11 using blue markers.
The bottom right figure of Fig. 11 shows the results obtained by Taubin fitting approach to estimate the circle parameters (Moon center and radius). The maximum distance from the estimated circle (worst residual) of the selected points was 2.9 pixels. The Taubin estimation gives a Moon center located at row 185.56 and column 187.34 and radius 218.03. Using this estimate, the set of points shown by red and blue markers in Fig. 13 over of Moon edge have been selected. These points are then used by CSF-LS to perform the nonlinear least-square estimation of Moon center and radius using circular sigmoid function. The new estimate implies a variation of −0.7, −0.5, and 0.4 pixels for row and column center coordinates and for radius, respectively.

ERROR ANALYSIS
This section is dedicated to the error analysis associated to the following sources: 1. Estimated position error. This analysis quantifies how the error in the measured Moon radius (in pixel) affects the Spacecraft-to-Moon distance computation (in km).
2. Estimated attitude error. This analysis quantifies how the error in the camera attitude knowledge (angle between quaternion) affects the Spacecraft-to-Moon direction (unit-vector).
3. Image processing errors. This Montecarlo analysis shows the distribution of the Moon radius and center estimated by the image processing software when the image is perturbed by Gaussian noise. This is done in the three cases of crescent, gibbous, and full Moon illumination.

Position and attitude uncertainty propagation
This section investigates the error in the Moon estimated radius because of the position estimation error. That is, how the Spacecraft position estimate affects the measurement of the Moon radius in the imager. It also investigates the error in the Moon center direction because of the attitude estimation error. That is, since the Spacecraft position estimate requires the evaluation of the Spacecraft-to-Moon vector in J2000 while it is measured in the camera reference frame, how the camera to J2000 change of coordinate affect the Spacecraft-to-Moon direction? Let us consider the Moon radius (in pixel) uncertainty first. Moon radius is directly related to the Spacecraft-to-Moon distance. The Spacecraft-to-Moon vector is r om = r m − r where r m and r are the Moon and Spacecraft position vectors in J2000, respectively. Vector r m is simply derived from time while the Spacecraft position vector, r, it can be considered a random vector variable with meanr and standard deviation σ r , that is, r ∼ N (r, σ rû ), whereû is a random unit-vector variable uniformly distributed in the unit-radius sphere. Neglecting the ephemeris errors, the uncertainty in the Spacecraft vector position coincides with the uncertainty in the Spacecraft-to-Moon distance, σ D = σ r .
The Moon radius in pixels in the imager is provided by where d is the pixel dimension (in the same unit of the focal length, f ), R m = 1, 737.5 km is the Moon radius, and D is the Spacecraft-to-Moon distance. This equation allows us to derive the uncertainty in the Moon radius in the imager The inverse of Eq. (41) provides an important information about the accuracy required by the radius estimation. This equation tells you how the radius estimate accuracy is affecting the Spacecraft-to-Moon vector length, and, consequently, how this estimate is affecting the estimation of the Spacecraft position. The Moon direction uncertainty (in pixel) is caused by the uncertainty in the attitude knowledge. Consider the attitude q ∼ N (q, 2σ qû ) with meanq and standard deviation 2σ q . Now the displacement with respect to the camera optical axis is c = f d tan ϑ, where ϑ is the angle between estimated Moon center and camer optical axis. Therefore, the radial uncertainty in the Moon center direction is provided by Image processing error Image processing error quantification requires performing statistics using a substantial set of images with known truth data. Since no such ideal database is currently available, the analysis of the resulting estimated parameters (Moon radius and center) dispersion is performed by adding Gaussian noise to the same image. The Moon illumination scenario plays here a key role since when the Moon is fully illuminated the number of data (pixels) used to perform the estimation will be twice as much as in the partial illumination case (under the same other conditions: distance to the Moon and selected camera). For this reason two scenarios have been considered: 1. Partially illuminated Moon (results provided in Fig. 21) 2. Full illuminated Moon (results provided in Fig. 22) The results of the tests done using partial illumination Moon are shown in Fig. 21. The image selected is real and it was taken on March 6, 2013 at 05:24:30 using a focal length f = 100 mm. The statistics of the results obtained in 1,000 tests by adding Gaussian noise with zero-mean value and standard deviation σ = 10 gray tone are given. The code estimates the Moon radius and center (row and column) by least-squares using circular sigmoid function (CSF-LS). The results show the three parameters estimated as unbiased and with a maximum error of the order of 0.2 pixels (with respect to the original picture). The histogram of the number of iterations required by the least-squares to converge shows an average of 8 iterations. The central-upper plot in Fig. 21 shows the ellipsoid distribution of the Moon center error due to the fact that the estimation uses just data on one side of the Moon, on the illuminated part. The results of the tests done using full illuminated Moon are shown in Fig. 22. Also in this case the image is real and was taken in Houston on February 25, 2013 at 22:33:00 CDT using a focal length f = 105 mm. The statistics of the results obtained in 1,000 tests by adding Gaussian noise with zero-mean value and standard deviation σ = 10 gray tone are given. The code estimates the Moon radius and center (row and column) by least-squares using circular sigmoid function (CSF-LS). The results show the three parameters estimated as unbiased and with a maximum error of the order of 0.01 pixels, one order magnitude more accurate than what obtained in the partial illuminated case. The histogram of the number of iterations required by the least-squares to converge shows an average of 3 iterations. The central-upper plot in Fig. 21 shows the unbiased Gaussian distribution of the Moon center error. These more accurate results are due to the fact that the estimation uses data all around the Moon edge.

CONCLUSIONS
This paper presents an image processing algorithm capable of determining the apparent location of the Earth or moon and its apparent diameter. This information can then be used for autonomous onboard trajectory determination in cislunar space. Tow method of removing potential outliers are introduced in order to increase the algorithm robustness. Preliminary analysis of the algorithm's performance is studied via processing both real and synthetic lunar images.