Lunar Crater Identification in Digital Images

It is often necessary to identify a pattern of observed craters in a single image of the lunar surface and without any prior knowledge of the camera’s location. This so-called “lost-in-space” crater identification problem is common in both crater-based terrain relative navigation (TRN) and in automatic registration of scientific imagery. Past work on crater identification has largely been based on heuristic schemes, with poor performance outside of a narrowly defined operating regime (e.g., nadir pointing images, small search areas). This work provides the first mathematically rigorous treatment of the general crater identification problem. It is shown when it is (and when it is not) possible to recognize a pattern of elliptical crater rims in an image formed by perspective projection. For the cases when it is possible to recognize a pattern, descriptors are developed using invariant theory that provably capture all of the viewpoint invariant information. These descriptors may be pre-computed for known crater patterns and placed in a searchable index for fast recognition. New techniques are also developed for computing pose from crater rim observations and for evaluating crater rim correspondences. These techniques are demonstrated on both synthetic and real images.


Introduction
Future lunar exploration missions are expected to rely on optical measurements (e.g., images from a camera) to navigate independently of Earth-based operators. Although the use of images for spacecraft navigation-called optical navigation (OPNAV)-is a well-established practice for conventionally navigated spacecraft [108,109], autonomous on-board OPNAV remains an emerging technology. Recent decades have witnessed great technological advancement and expanding acceptance of autonomous vision-based navigation for the exploration of other celestial bodies (e.g., Moon, Mars, asteroids). These advancements in real-time, on-board OPNAV are exemplified by the technological progression from simple estimation of lander velocity with the Mars Exploration Rover's DIMES system in 2004 [19] to autonomous feature tracking that will soon be demonstrated on the OSIRIS-REx mission to asteroid Bennu [105] and during landing of the Mars Perseverance (Mars 2020) rover [64].
This work focuses on OPNAV for lunar missions using digital images captured by a conventional camera (images from a pushbroom camera [47] are a different problem and are not discussed here). Images from a conventional camera follow the geometry of perspective projection. How to best use such images of the Moon for navigation depends on many factors-with distance and lighting often being the two most important considerations. It is usually best to use horizon-based OPNAV [21,26,108] when far away from the Moon and when a large portion of the lunar disk is contained within the image. As the vehicle gets closer to the Moon (e.g., low lunar orbit, lunar lander descent), it becomes more appropriate to observe specific landmarks on the lunar surface. These landmarks could be specific types of morphological features (e.g., craters) or a patch of unique-looking terrain.
The principal difficulty with landmark-based OPNAV is the need to match points in an image to points in an onboard map. This has led to the continued use of horizon-based OPNAV at close distances where it would otherwise be better to use surface features. It has also led to navigation with unknown landmarks [8,23,80,137] or visual odometry [27]. Unfortunately, horizon-based OPNAV is not always practical (especially at close ranges) and unknown landmarks do not generally provide full state observability. Consequently, the ability to match surface features to a map is a required capability for autonomous vision-based navigation for many missions.
Many landmark-based OPNAV algorithms rely on the real-time rendering of an onboard digital elevation map (DEM) [3,44,64,105]. The standard approach is to render the expected appearance of landmark patches and then compare this to regions of the navigation image, usually by means of a 2D cross-correlation. This, however, requires the vehicle to carry onboard 3D models for each landmark patch (and, sometimes, complete DEMs for certain portions of the body) and also have the ability to render these patches in real-time. Further, the rendering step requires a priori state knowledge, thus precluding this technique from supporting lost-inspace landmark-based OPNAV. More troubling, however, is that the rendering-based template matching approach feeds navigation state data into the measurement generation process, creating scenarios where corrupted navigation states render otherwise good images ineffectual-potentially leading to unnecessary filter reinitializations, trajectory aborts (e.g., during lunar descent), or other undesirable events.
Landmark-based OPNAV using craters [18,49,86,110] is an alternative strategy that, when properly implemented, does not necessarily require a priori state data. Furthermore, OPNAV with craters does not require onboard DEMs or an onboard rendering capability. Instead, such systems require an image processing algorithm to detect craters and a pattern matching algorithm to match these observations to a catalog of known craters. This manuscript focuses on the second part (crater pattern matching) of this problem.
There are two general classes of crater matching algorithms: (1) tracking and (2) lost-in-space. In the case of crater tracking, it is possible to use a priori state information to predict where catalog craters should appear in an image. These predictions are used to associate observed craters with model craters through a variety of schemes of varying sophistication. We generally find that explicit crater matching is undesirable and more robust tracking performance is achievable using probabilistic techniques (as exemplified by anonymous feature tracking (AFT) [88]). In the case of lost-in-space crater matching, no a priori state information is available. Lost-in-space algorithms, therefore, must be able to recognize a crater pattern from its appearance in an image regardless of camera pose (i.e., camera position and attitude). Such a capability is necessary for filter (re)initialization. It is also useful in providing state-independent measurements (if required) or a state-independent verification of crater tracking correspondences. This work focuses on solving the lost-in-space crater identification problem.
We briefly note that, while autonomous spacecraft navigation was the original motivation for this work, crater pattern recognition is equally useful for the registration of scientific images [138].
In this manuscript, we provide a comprehensive treatment of lunar crater pattern recognition in digital images formed by a conventional camera (i.e., under perspective projection). We begin by outlining the philosophical framework for such a system (Section 2) and then review important background material (Section 3). This is followed by the detailed mathematical development in Sections 4-9. Numerical results and examples are shown in Section 10. This work contains a number of key results, which are summarized here: 1. Almost all lunar impact craters are elliptical, and many are nearly circular (Section 3.1). Rather than being a heuristic design choice for easy crater identification, we discuss why craters are necessarily elliptical in shape from a perspective of planetary science and impact mechanics. Recent advancements in our knowledge of the lunar crater population suggests that a circular crater assumption (common amongst crater identification algorithms) is not well-supported by the data for small craters (diameter less than about 30 km). 2. Invariant theory is the proper mathematical framework for recognizing a pattern of crater rims in an image (Section 5). Elliptical craters on the lunar surface project to elliptical features in an image (Section 4). For some patterns of craters, there exist algebraic quantities that remain unchanged (i.e., invariant) regardless of the camera pose and that may be computed from just the projected crater rims in a single image. Such quantities we call projective invariants.
Invariant theory provides the means by which we may determine the existence (or not) of these projective invariants. 3. There are no projective invariants for arbitrarily placed conics (e.g., elliptical crater rims) on the surface of an arbitrarily shaped body. We provide the first known proof of this fact (Section 5.4.3). This result is significant, as it produces a theoretical obstacle for solving the lost-in-space crater identification problem about very irregular bodies, such as some asteroids and comets. Fortunately, the Moon is not an arbitrarily shaped body. 4. Projective invariants exist for arbitrarily placed conics (e.g., elliptical crater rims) lying on a nondegenerate quadric surface (e.g., sphere, ellipsoid), which is an excellent approximation for the shape of the Moon. A different set of invariants exists for conics lying on a common plane. The existence of invariants for conics on a nondegenerate quadric surface is a novel result (Section 5.5), while invariants for conics on a plane is a known result (Section 5.2). Since the Moon is nearly spherical on the regional/global level and nearly planar on the local level, projective invariants generally exist for craters lying on the surface of the Moon. In both cases (global and local), we show how to efficiently compute these invariants from only the projected conics (i.e., the contour of the crater rim) that appear in a digital image (Section 6). 5. There are no algebraically independent invariants for conics (e.g., elliptical crater rims) on a nondegenerate quadric surface or on a plane beyond the ones developed in this work. All other invariants a researcher may conceive for one of these two cases may be a written as an algebraic function of the invariants discussed herein. Thus, our invariants describe all the independent information about a crater rim pattern that is pose invariant and useful in constructing a descriptor that may be indexed. We provide the first known proof of this fact for the case of conics on a nondegenerate quadric surface (Section 6.1.3). Given these findings, we suggest future investigation of invariants for patterns of lunar crater rims be focused on ease of computation or numerical stability rather than attempting to extract additional independent information from the crater rim contours (since there is no additional independent information to be had that remains unchanged with camera viewpoint). 6. The projective invariants for a triad of lunar craters may be used to form a feature descriptor that is insensitive to camera viewpoint (Section 7). These descriptors may be computed for known crater patterns and stored in a searchable catalog (i.e., an index). Matches to an observed crater pattern in any image are found by simply finding the nearest neighbor for its descriptor in the index. In cases where it is desirable for the matching to be permutation sensitive, the descriptor is nothing more than the projective invariants concatenated into a vector. Alternatively, in cases where we wish the matching to be insensitive to crater ordering, we apply invariant theory a second time to construct projective and permutation (p 2 ) invariants. This is a novel result for crater pattern descriptors. 7. Lunar crater identification is a multi-scale pattern recognition problem. Local, regional, and global crater patterns are built from individual craters of different (increasing) diameters and occur over different (increasing) extent on the lunar surface. For example, a 1 km crater may contribute to a local crater pattern but not a global crater pattern, while a 100 km crater may contribute to a global crater pattern but not a local crater pattern. To address this challenge, we present the first known use of Hierarchical Equal Area isoLatitude Pixelization (HEALPix) to construct a hierarchy of crater catalogs (Section 8). 8. A substantial amount of navigation information is contained within the shapes of projected crater rims in an image. Many existing crater identification algorithms, however, compute camera pose using only the center coordinates of the craters-thus ignoring a great deal of actionable information. We discuss how to compute the camera position using the contours of projected image conics, rather than just the conic center coordinates. Our solution is non-iterative and finds the camera position in the least squares sense when observing d ≥ 2 craters (Section 9.1). 9. Crater match hypotheses are verified by comparing the observed crater rims to what is expected were the hypothesized match true. Most existing crater identification algorithms perform this verification by comparing only the coordinates of crater centers, which generally requires ≥ 5 craters locations to agree. Instead, we develop a novel distance metric as a measure of how dissimilar two ellipse contours are from one another (Section 9.2) and use this to compare the observed and expected crater rims. Our metrics satisfy the three classical axioms for a metric (minimality, symmetry, triangle inequality) along with a fourth requirement of similarity invariance. Using the entire crater rim permits pattern verification with just the three craters used to perform the index match and does not require additional craters for verification. Thus, as compared to past methods, the new distance metric permits crater identification over more sparsely cratered regions of the Moon. 10. There is no single crater identification algorithm that is always best. The correct solution depends on orbital regime (or descent trajectory), camera specifications, on-board computational and memory resources, operational cadence, and other mission-specific parameters. There are critical crater identification design decisions to be made for the crater catalog (raw source of data), index scale (local, regional, global, or a combination), index organization, and invariant descriptor type (projective invariants or p 2 invariants). Therefore, rather than present a single crater identification algorithm, we suggest a framework solving this problem and provide the reader with tools for each step. Once understood, the elements in this manuscript may be used to identify craters for a wide diversity of mission types. A few illustrative examples are provided to enhance understanding of how this might work.
Finally, we note that this work relies on mathematics and techniques that may be obscure to the average astrodynamicist and spacecraft navigator. Therefore, we adopt a rather explanatory approach to facilitate understanding, dispel common misconceptions, and encourage adoption of the techniques herein.

A Framework for Lost-in-Space Terrain Relative Navigation (TRN)
The principal difficulty with recognizing 3D objects from their projection in a digital image is that the appearance of these objects changes with camera viewpoint and illumination geometry. The recognition of landmarks on the lunar surface inherits this generic computer vision difficulty. In many cases (including lunar landmark identification), we may circumvent this problem entirely by choosing to represent objects with a descriptor that remains the same (is invariant) regardless of camera viewpoint. Such an object recognition philosophy was influentially advocated for in [97] and [147], and we adopt this philosophy here.
There are a variety of object attributes that may remain unchanged with camera viewpoint and illumination geometry, including color, texture, and geometric structure. Finding these invariant attributes, however, requires a great deal of care. Of note, we know that no such invariant exists for an arbitrary patch of terrain with constant albedo and with an approximately Lambertian reflectance [17]. This suggests that terrain patches (or maplets)-as is often used in correlation-based TRN [3,44,64,105]-may not be well-suited for lost-in-space TRN, though this is an important topic of future work. We also observe that popular image feature descriptors-such as SIFT [83], SURF [7], BRIEF [78], ORB [124], and others-do not describe the type of invariance required for lost-in-space TRN. They are, by construction, not the appropriate tool for the problem considered in this manuscript. There are two primary reasons for this. First, most of these descriptors are only formally invariant for a similarity transformation (translation, rotation, scaling), though many have been shown to be robust for a modest amount of affine shear. While there do exist feature descriptors that are formally affine invariant [92], it was convincingly argued by Lowe [83] that such invariance is often undesirable within the context of feature descriptors. Regardless, none of these descriptors are especially robust to extreme projective transformations-making them more appropriate for matching features between two similar images than for matching a single arbitrary image to a prebuilt database of features on a 3D object. The second (and more significant) problem with classical feature descriptors for lost-in-space TRN is that none of them are illumination geometry invariant, as discussed at length in [27]. In brief, all of these feature descriptors work by finding and describing regions of unique 2D patterns in an image. For TRN above an airless body (e.g., Moon, asteroid), the 2D pattern observed in an image is formed by the way that light is reflected off the terrain and towards the camera. When the lighting geometry changes, the 2D pattern in the image changes-thus the feature descriptor for the same surface patch will change. Hence, the feature descriptors are not illumination invariant.
An alternative to landmark patches and feature descriptors is to look for explicit geometric features. Crater rims are one such geometric feature that can be easily recognized and consistently localized, even when observed from different camera viewpoints and under different lighting conditions. We find the geometry describing the pattern of multiple crater rims to be one of the most accessible attributes from which viewpoint invariant information may be constructed.
For any given object geometry, such as a pattern of crater rims, there may be a variety of viewpoint-invariant properties-but not all of these necessarily provide equal perceptual salience [63]. Thus, at a minimum, we seek non-constant invariants 1 that can effectively differentiate between many different crater rim patterns. Such non-trivial invariants do not always exist for craters lying on the surface of arbitrary 3D celestial bodies, though they do exist for a body such as the Moon. Even when invariants do exist, finding a way to compute them is not always straightforward. However, as we show in this manuscript, the additional effort to develop these invariants produces powerful pattern recognition capabilities.
The identification of perceptually salient invariants under the action of perspective projection (e.g., a digital image of the Moon collected with a conventional camera) allows us to quickly distinguish between two different objects (e.g., two different crater patterns) possessing two different invariant descriptors. Therefore, if we compute the invariants for a particular lunar crater pattern, these invariants may be concatenated into a descriptor for that pattern and stored in a catalog (or index) for future comparison. We repeat this procedure for each known crater pattern preflight and produce a very large index of patterns and their corresponding descriptors. Then, when an image is acquired in-flight, the invariants may be computed for an observed crater pattern and concatenated into the same type of pattern descriptor that we stored within the index. Therefore, the index may be queried for known patterns with descriptor values similar to that of the observed crater pattern. Efficient data structures allow such queries to be performed very fast, usually in O(log n) time.
Matches from the index may be used to construct crater match hypotheses, which must be verified before acceptance. The primary means for verification is to use a hypothesized crater correspondence to compute the camera's location. This hypothesized camera location is used to reproject expected crater rims into the image, which are compared to the observed crater rims. If the reprojected pattern matches the observations, the match hypothesis is accepted. If not, the match hypothesis is rejected, and we attempt a different hypothesis.
This framework is summarized in Fig. 1. The remainder of this manuscript is dedicated to providing the details for each step of this process. Fig. 1 The crater identification framework has two major components: index construction and pattern recognition. The index construction step (left) is often time-intensive since we must loop through all possible triads and then store them in an efficient data structure. The crater patter recognition step (right) is fast and is performed for each image

Observations on Lunar Crater Morphology and Size-Frequency Distribution
Lunar crater identification algorithms often begin with the assumption that the crater rims are either circular or elliptical in shape. It is also generally assumed that craters are well distributed about the Moon and that the crater catalog is static (new craters are not being added or removed). We find, however, that many recent crater identification algorithms are built on assumptions not supported by a modern understanding of the lunar crater population-and, therefore, are not likely to perform well across the full diversity of orbital regimes. A thorough understanding of lunar crater morphology and size-frequency distributions must precede the development of a computational algorithm to recognize patterns of these features.

Crater Catalogs
To recognize crater patterns in images of the lunar surface, we must first know where the craters are located and how they are shaped. Such data is contained in lunar crater catalogs. In most cases, however, these catalogs were built for planetary science purposes and not for image registration purposes-with the specific scientific question motivating the catalog construction often influencing the contents of the catalog. This must be considered before repurposing a catalog for spacecraft navigation, since not all lunar crater catalogs are appropriate for building an index for crater pattern identification.
There are a variety of lunar crater catalogs available; e.g., [52,112,119], and others [126,140]. These catalogs were created with differing scientific objectives and from different raw data sets-thus, there should be no expectation of a one-to-one correspondence between their entries. For example, the minimum crater diameter is different for each catalog (and some catalogs even focus on very specific size ranges). As another example, some catalogs include all crater-like objects (e.g., [119]), while others only include craters the authors are confident originated from an impact event (e.g., [112]). Moreover, some catalogs attempt to only count the craters produced by direct impacts, thus ignoring those craters thought to be produced by secondary impacts (i.e., ejecta from the primary impact). This highlights the importance of selecting the crater catalog that best matches a particular mission's crater identification needs.
Crater catalogs are a useful tool for studying the statistical properties of the entire lunar crater population. To obtain global and comprehensive coverage, we usually exchange the detailed understanding of individual craters for a general understanding of all craters. Clearly, a detailed study of any single crater would yield a more sophisticated understanding of that crater-and likely produce crater information (e.g., size, shape, depth) that differs slightly from its corresponding entry in a global crater catalog.
The most comprehensive crater catalog to date is that of [119], containing about 1.3 million lunar impact craters having a diameter larger than about 1-2 km. This is the database we choose to use in this work.

The Elliptical Crater Assumption
The natural shape of an impact crater is an ellipse. Almost all lunar craters are nearly elliptical, and many craters are nearly circular [14]. The same is true for crater populations on other large celestial bodies [54]. Clearly, no crater rim is a perfect ellipse or perfect circle and this is only an approximation of the general shape.
The are a variety of factors that determine the size and shape of an impact crater. With regards to the shape of lunar crater rims, the most important factors are the impactor velocity, impactor size, and lunar surface properties. Lower impact angles, φ (see Fig. 2), tend to produce craters of higher ellipticity-where we define ellipticity, = a/b ≥ 1, as the ratio of crater rim semi-major to semi-minor axis. The sensitivity of the crater ellipticity to this impact angle is a function of the cratering efficiency (the ratio of crater diameter to impactor diameter) [14,39], which acts as a model surrogate for impactor size and lunar surface properties. An empirical fit to hydrocode simulations and laboratory experiments shows excellent agreement between the impact angle at which = 1.1 occurs (φ =1.1 ) and a simple power law relation [39], where D 90 is the crater diameter for a vertical impact and L is the impactor diameter (i.e., D 90 /L is the cratering efficiency for a vertical impact). Thus, an impact with φ < φ =1.1 will produce a crater with > 1.1. The cutoff angle φ =1.1 is generally between 3-5 deg (for sand, D 90 /L ∼ 60) and 30 deg (for brittle rock, D 90 /L ∼ 4) [29,39,91]. Given the rather oblique impact angle required for an elliptical crater, it follows that a large percentage of lunar craters should be nearly circular. For an airless body such as the Moon, the cumulative distribution function (CDF) for the impact angle when the impactors are assumed to arrive from random directions is [130] p(φ ≤ φ 0 ) = sin 2 φ 0 (2) This model suggests that about 7.6% of the crater population should have an impact angle less than 5 deg and that about 25% of the crater population should have an impact angle less than 30 deg. The database of [119], however, found a higher percentage of elliptical craters than might otherwise be expected and proposed a few hypotheses for this observation. Of particular note is that Robbins's database includes both primary and secondary impact craters, whereas many ellipticity studies (e.g., [14]) only consider primary craters. Regardless of the underlying physical cause, the newly available global distribution of crater ellipticity from [119] makes clear that a circular crater assumption is not valid at the local level (i.e., when viewing craters of diameter smaller than about 30 km). This finding is significant, as many recent crater identification algorithms assume a circular crater model-a choice that is not well supported by the data. To see this, consider the ellipticity distribution as a function of crater size as shown in Fig. 3. We immediately see that over half of the entire crater population has an ellipticity greater than 1.1 (the threshold planetary scientists usually use to define an "elliptical" crater [14,39]). Many of these craters, however, are not morphologically well preserved-and, therefore, are not good candidates for crater-based navigation since they will not produce good ellipse fits. If we restrict ourselves to well-defined craters whose catalog shape is supported by over 90% of its circumference, we see a substantially larger percentage of the craters are nearly circular, especially when considering craters larger than 30 km.

Fig. 3
Lunar craters with well-defined rims tend to be less elliptical than the general population of craters. Plots show the fraction of craters above a specified ellipticity for every catalog crater (top) and for only those with well-defined crater rims (bottom). A visualization of crater ellipticity is shown above the plots for intuition and context. Results are based on post-processing of catalog data from [119] Topography or layered strength differences can create other irregular crater shapes [89], with a few examples shown in Fig. 4. On the Moon, simple craters ( 10-20 km in diameter) often have bowl-shaped morphologies, with small flat floors and few internal topographic features. Some may have deposits on the floor that originate from mass wasting or ponding of impact melt. The rims of simple craters are generally uplifted (unless significant degradation has occurred, as in the case of old simple craters), and the crater is surrounded by an ejecta blanket. More energetic impacts create larger, more complex craters. Complex craters ( 10-20 km) have more varied morphologies, including central peaks, terraced walls, and flattened floors. Central peaks are formed by rebounding of material from deeply buried rocks, making these features especially interesting for studying materials excavated from lower in the crust. Larger complex craters can often have more than one central peak structure. Surrounding rocks in the wall of complex craters can also collapse, forming terracelike structures. As impact features increase in size, they also increase in complexity. Craters larger than 300 km are usually classified as basins.

Lunar Crater Size-Frequency Distribution
Analyzing crater populations plays a major role in establishing planetary chronologies and in dating geologic units. There is general consensus that the lunar cratering projectile flux was relatively constant within the last 3 Gyr [102,133], with the period before 3 Gyr experiencing variable impact rates as a result of the hypothesized late heavy bombardment.
Cumulative crater size-frequency distributions (SFDs), which describe the number of craters greater than a given diameter per measured area, are the usual data product used to estimate the age of planetary surfaces. Crater SFDs are conventionally determined by binning measured craters based on their diameters [46,55], though modern practitioners are transitioning to more formal statistical approaches (e.g., using empirical density functions) [120]. These SFDs are often used to produce production functions that describe the rate at which new craters occur.
The most commonly used production function for estimating lunar surface ages is the Neukum Production Function (NPF) [101,102]. The NPF models the formation rate of craters between 10 m and 300 km, and is generally used to estimate ages for geologic units on the Moon [102]. The NPF is shown in Fig. 5 and is given by where D is the crater diameter, N is the number of craters of size D or smaller per unit area per unit time, and the coefficients a k (for m = 11) are given in [102]. The NPF suggest that that we should expect the production of new craters on the lunar surface to occur at a rate of 1.7×10 −5 new craters larger than 4 km per km 2 /Gyr. For the Moon, this results in approximately 635 new craters per Gy (or about 1.7 × 10 −9 craters/day). Assuming impacts follow a Poisson distribution, we may straightforwardly compute the probability of an impact creating a crater larger than 4 km during a specified period of time. If a crater catalog must be valid for a period of days to a few years, the probability of a new crater appearing is essentially zero. Thus, it is quite reasonable to assume a static crater catalog.
We note that recent work analyzing new impact craters as detected by the Lunar Reconnaissance Orbiter (LRO) Camera [133] found that the current cratering rate may potentially be higher than what the NPF predicts. Regardless, the expected discrepancies are small within the present context and the static crater catalog assumption remains valid.

Crater Detection
Successful crater detection is a prerequisite for crater identification. While not the focus of the present work, a brief review of crater detection algorithms (CDAs) is appropriate. There have been a number of recent surveys of CDAs for both planetary science [33,125] and spacecraft navigation [142] applications. We observe here that the crater detection requirements for planetary science and navigation are often quite different-meaning algorithms developed for one application are not necessarily transferable to another application. There are two primary differences. First, from a scientific standpoint, the CDA objective is often to find all of the craters. For navigation, however, the objective is not to find every crater, but only enough to navigate. Any crater found that is not in the database becomes clutter that must be rejected. Thus, an algorithm that returns every crater is not necessarily better for navigation, and we are more interested in algorithms that preferentially return only well-modeled craters that exist in our on-board catalog. The second primary difference between scientific and navigation CDAs is suitability for real-time use. CDAs for autonomous navigation must be fast enough to run in real-time on a space-qualified computing platform and robust enough to not require human supervision. No such requirements are generally placed on CDAs for scientific analyses.
With well over 100 different algorithms contained in the CDA surveys [33,125,142], we find no need to discuss the topic further here. The interested reader is directed to these surveys and the references therein.

Crater Identification
Lost-in-space crater identification from a single image is a pattern recognition problem. When challenged with an image containing a set of observed craters, the task is to find the corresponding craters from within a large catalog of known craters.
This problem is not new, and numerous attempts have been made in the past to recognize crater patterns. However, the lack of an existing invariant for non-coplanar ellipses has forced past work to make numerous simplifying assumptions or employ ad hoc recognition schemes. These restrictive assumptions are relaxed by the developments introduced in this work. We now review some of the most notable prior work in crater identification and discuss their drawbacks using insights from invariant theory. As the reader will see, miscounting of independent invariants is a pervasive problem in the crater identification literature-leading to descriptors that are incomplete, redundant (i.e., not all elements are independent), or both. This recurring problem is one of the principal motivations of the present work.
Hanak [49,50] uses an adaptation of star identification algorithms [95,127] to recognize crater triads-an interpretation of the crater identification problem that has become influential within the spacecraft navigation community. Hanak chooses to index the crater pattern using the triangle's interior angles, which is not an invariant under general perspective projection. Therefore, the approach of [49,50] operationally restricts usability to a nadir pointing camera with a modest field-of-view, which restricts the crater pattern to a similarity transformation where inter-crater angles remain invariant. If the camera is not nadir pointing, or if the spacecraft is far away from the Moon and craters are not nearly coplanar, the descriptor is no longer pose invariant and matching performance may be poor. Once potential matches are made, the authors of [49] check these hypotheses against four additional metrics (three ratios of crater diameters to inter-crater distances, and a ±1 flag indicating if the triangle is clockwise or counterclockwise going from shortest to longest leg), which are also not projective invariants. We observe, therefore, that the authors of [49,50] develop a six element descriptor, with none of the elements being a projective invariant. We contrast this to the seven known invariants for three coplanar conics. It is possible, therefore, to build a seven-element descriptor (see Section 6.2.2) that holds more descriptive power than [49] and that is also a projective invariant.
Yu, et al., assume coplanar craters under affine transformation [145], and develop an invariant descriptor based on the ratio of crater areas. The ratio of areas of two closed coplanar curves as invariant for affine transformations is discussed in [41]. Identification using only area ratios, however, neglects the spatial relationships between the craters and significantly reduces the descriptiveness of these invariants. Like the work of Hanak [49], the methods from [145] are also limited to nearly nadir pointing images of nearly coplanar craters.
Cheng, et al., develop a crater identification scheme by matching pairs of coplanar craters [18,20] using invariants from [41] and [115]. A pair of coplanar conics possess two unique invariants, which may be computed preflight and indexed into a catalog for in-flight matching. The descriptive power of only two conics, however, is low-so multiple pairs must be compared, which is straightforward to accomplish. Additionally, the constraint of using coplanar conics limits the flight regimes where this technique works well. Despite these minor drawbacks, we consider the work of [18,20] to be the most theoretically informed lost-in-space crater identification approach developed to date.
Park, et al., consider a pattern of three coplanar craters [110] by looking at invariants of six coplanar points, where the six points are formed by the intersection of the triangle connecting the crater centers with the crater conics. These six coplanar points are split into six possible groupings of five points, from which the two unique invariants for each set of points may be computed using results from [41]. Then, for each set of points, these two invariants are made insensitive to point ordering using the projective and permutation invariants (also called p 2 invariants) from [76]. The pair of invariants for five coplanar points from [41] produce five p 2 invariants [76]. The authors of [110] concatenate the five p 2 invariants for each of the six sets of points to create a 30-element descriptor (5×6 = 30) for the crater triad. However, as discussed in [76], only two of the five p 2 invariants for a set of five coplanar points are independent (which should be expected, since making the two invariants from [41] insensitive to point ordering can not add new invariants). Therefore, the six sets of five points are really described by vector of 6 × 2 = 12 invariants, such that the 30-element descriptor from [110] has 18 redundant elements that do not enhance discriminative power. Moreover, the six sets of five points are not constructed from unique points (they are the six possible five-point combinations of the same six points), thus the remaining 12 invariants are not all independent. We observe that six coplanar image points possess only four independent projective invariants. We also observe that d ≥ 2 coplanar conics possess only 5d −8 independent projective invariants (seven invariants for a triad of craters), as discussed in Section 5.2. Consequently, the 30-element descriptor of [110] has only four independent elements, 26 redundant (dependent) elements, and is missing three independent invariants. We note here that the seven invariants for three coplanar craters may be computed directly from the observed image conics as discussed in Section 6.2, without the need to compute sets of points. Likewise, invariants for clusters of four or more coplanar craters may also be computed directly [53]. The practical mechanics of computing such invariants are discussed in Section 6.2.2. If one wishes to build a permutation invariant descriptor (e.g., with p 2 invariants), we discuss how to do this for triads of craters in Section 7.
Maass, et al., propose an assortment of crater identification techniques [86], though we only discuss their lost-in-space technique here. They assume circular craters, thereby allowing explicit estimation of the crater normal in the camera frame (there are two possibilities [22,129]). Given a set of three craters and their normals, the authors of [86] consider the 2 3 = 8 possible configurations and choose the one where the normals are most similar (which presupposes a nearly planar configuration). More craters can be added to the pattern by iteratively solving a binary global optimization problem, with the cost function being the curvature energy of a cubic spline through a Delaunay triangulation of the crater centers. This resolves the 3D crater configuration (or reconstruction) to an unknown scale. Thus, Maass, et al., define a pattern descriptor using a pair of independent interior angles of a 3D crater triad, which we observe is a projective invariant since it is a reconstruction of the 3D geometry (and not the triangle in the image). Maass, et al., augment this check with the two independent ratios of crater radii. Pattern matches are verified by comparing the reprojected centers of additional craters. This usually requires two additional craters, thus requiring an image to possess five craters to yield a successful match. The method of [86] just described has three major drawbacks. First, we know from Section 3.1.2 that the circular crater assumption is often not valid, with 25%-50% of small craters having an ellipticity > 1.1. This effects the ability to accurately estimate the surface normal direction (e.g., a nadir-pointing image of an elliptical crater with = 1.1 would have a normal direction error over 20 deg if it was incorrectly assumed circular). Second, the two parameter pattern descriptor (or four parameters if you count the radii ratios, though these do not appear to be explicitly used in the search) are not maximally descriptive-we know there to be seven invariant parameters for a triad of nearly coplanar craters (see Sections 5.2 and 6.2.2). Thus, while the descriptor is indeed pose invariant (if the craters are circular), there are at least three independent pieces of information that remain for describing the crater pattern. Third, by verifying matches using only crater centers, more craters are required. Verifying using the crater rim reprojection (see Section 9.2) reduces the number of craters required.
There are a variety of other crater matching schemes [31,68,84], all of which make assumptions similar to the above methods.
This work is most similar to that of [18,20], though we relax the requirement that craters be coplanar-thus enabling a substantial expansion of the orbital regimes where crater matching may be performed. We also allow for the simultaneous consideration of more craters by using triads instead of pairs.

Mathematical Representation of a Crater
Since impact craters are known to be elliptical in shape, many lunar crater databases store the ellipse parameters for the best-fit ellipse to the rim of each cataloged crater. The most common database parameterization is the (1) latitude/longitude of the crater center {ϕ, λ}, (2) lengths of crater ellipse semimajor and semiminor axes {a, b}, and the orientation of the crater ellipse ψ (often measured counterclockwise from East). This is the parameterization used in the Robbins database [119].
In this work, we assume that a crater is a planar feature defined by the 2D elliptical curve that best describes the crater's rim. Circular craters are a special case (a = b) that is naturally handled without special treatment.
Let ξ M be a three-dimensional (3D) point in the selenographic (i.e., Mooncentered, Moon-fixed) frame, which may be represented in homogeneous coordinates as A crater is a planar feature, so let the plane P i containing crater i be described by the 4 × 1 vector π i , such that Here, P n denotes projective space, which is discussed in greater detail in Ref. [51] and Ref. [43].

Two-dimensional (2D) Crater Description
Within the plane P i , define a 2D coordinate system with origin at the ellipse center and aligned with the ellipse principal axes (see Fig. 6). In this case, we arrive at the canonical form of an ellipse, where the point [x , y ] lies on the ellipse if where a is the ellipse semi-major axis and b is the ellipse semi-minor axis. Now, define another 2D coordinate system within P i with a different origin and different coordinate axis orientation (see Fig. 6). Let the center of the ellipse be located at [x c , y c ] and the ellipse principal axes be rotated by an angle ψ. A point [x, y] in this coordinate system is related to a point [x , y ] by a Euclidian transformation, Substituting this result into Eq. 7 yields the implicit equation

Fig. 6 Visualization of ellipse parameterization in P i
Observe that Eq. 9 is the generic equation for a conic, which happens to be an ellipse when B 2 − 4AC < 0. A simple calculation with Eqs. 10-12 shows that B 2 − 4AC = −4a 2 b 2 , which will always satisfy the ellipse inequality constraint when a ≥ b > 0. The reverse mapping (implicit coefficients to explicit parameters) is straightforward and left to the reader (though the expression may be found in [25]). The implicit equation from Eq. 9 may be written as a quadratic form using homogeneous coordinates. Therefore, letting a 2D point in the crater plane P i be given by which may be written in homogeneous coordinates as we observe that Eq. 9 may be rewritten as The expression in Eq. 18 describes a conic locus, which is the locus of 2D points forming the path of the conic (see Fig. 7). The matrix C describing a conic is a 3 × 3, real-valued, symmetric matrix of ambiguous scale. The ambiguous scale means that C has only 5 degrees-of-freedom (which is also obvious from Eq. 9) and that C and kC (with k being a real-valued, non-zero scalar) describe the same conic. If C is full rank (as it is with an ellipse), we have a proper conic-otherwise it is degenerate Fig. 7 Visualization of simple conic locus (black) and lines belonging to the conic envelope (red) (either a pair of lines or a double line). Furthermore, if the conic is an ellipse, C is indefinite-always having two eigenvalues of one sign and a third eigenvalue of the opposite sign, with none of the eigenvalues being zero (since it is a proper conic).
Since C is full rank for an ellipse, det(C) = 0 and where C * is the adjugate matrix of C.
In projective geometry, points are dual to lines [51,128]. A 2D point written in homogeneous coordinates (x) lies on the line if it satisfies the constraint Tx = 0. Clearly,x and can be exchanged in this relation-hence the point-line duality. To apply this duality property to conics, first observe that the line ∝ Cx is tangent to the conic C ifx is a point on the conic [51]. For a proper conic, one may therefore computex ∝ C −1 . Substitution into Eq. 18 yields and, since C is symmetric and C −1 ∝ C * , where is a line tangent to the conic. The utility of using adjugate matrices instead of inverses will become readily apparent as we proceed. The expression in Eq. 22 describes a conic envelope, which is the family of tangent lines that encapsulate the ellipse (see Fig. 7).

Three-dimensional (3D) Crater Description
The elliptical crater rim is a 2D feature in 3D space. Thus, as we will show, the crater catalog data may be used to construct both the crater's plane and the disk quadric describing the 3D crater. The disk quadric is the natural and mathematically rigorous way to represent a 3D conic feature.
Consider catalog crater i with a geometric center at lunar latitude ϕ i and longitude λ i . The 3D selenographic position of the crater center is where ρ i is the distance of the plane P i from the center of the Moon. The crater's center point is assumed to lie in the crater plane, π T ip (c) M i = 0. The Moon is very nearly a sphere. Therefore, to an excellent approximation, the crater normal is parallel to the selenographic crater center location p (c) M i . If the celestial body were an oblate spheriod, the local "up" direction may be computed through Thus, defining the columns of T E i M according to we arrive at a relation between the 2D ellipse pointx E i in the crater plane [see Eq . 17] and its corresponding 3D point in the selenographic framē and where S is the 3 × 2 matrix Consider a quadratic surface in P 3 (e.g., a sphere, an ellipsoid). Such surfaces are generally called quadrics, with a 3D pointξ lying on the surface if where Q is a 4 × 4 symmetric matrix describing the surface. This is the 3D analog to the 2D expression in Eq. 18. The 3D surface defined by Eq. 36 is called a quadric locus. The 3D conic locus (a curve) cannot be represented by a quadric alone. Instead the conic defines a proper quadric cone (where Q is a 4 × 4 matrix of rank 3) and the 3D conic is formed by the intersection this quadric cone with the plane P i . Therefore, define the Moon-centered quadric cone (i.e., a cone with vertex at the center of the Moon) as such that the conic locus describing the elliptical crater is formed by the conic section Just as points and lines are dual in P 2 , points and planes are dual in P 3 (a fact that should be evident from Eq. 6). It follows therefore, that one may construct a dual quadric, Q * i , such that the quadric envelope is given by where π is a plane tangent to the quadric. The quadric envelope provides a more natural way to describe a 3D conic. The quadric envelope for a 3D conic is generally called the disk quadric [128], and defines all the planes tangent to the conic curve (see Fig. 8). This may be visualized as all the planes tangent to an ellipsoid x 2 /a 2 +y 2 /b 2 +z 2 /c 2 = 1 and letting c → 0, thus resulting in a 3D disk (i.e., the 3D ellipsoid is collapsed to a plane and resembles a pancake or dinner plate). The disk quadric is, therefore, defined by a 4 × 4 matrix Fig. 8 Visualisation of the disk quadric that defines all of planes tangent to D i lying in the plane P i . A few example tangent planes are illustrated in red Q * i that has rank 3. We may compute the disk quadric directly from the crater's conic envelope as which is a 4 × 4 symmetric matrix of rank 3. The disk quadric is clearly a more direct way to represent the 3D conic than the intersection of a proper quadric cone and a plane. Thus, the disk quadric is the representation of choice when one must concisely describe a 3D conic, such as a crater located on the surface of the Moon.
Each 3D conic (as defined by the disk quadric) corresponds to a symmetric 4 × 4 matrix of rank 3, up to a scalar. The symmetric 4 × 4 matrices form a 10-dimensional vector space, so a nonzero 4 × 4 symmetric matrix corresponds to a point in P 9 . The constraint that the rank is ≤ 3 defines a hypersurface in the projective space P 9 given by the vanishing of the determinant. The set Z of all symmetric 4 × 4 matrices of rank 3 is an open subset in this hypersurface that parametrizes the set of conics. The dimension of Z (and the hypersurface) is 9 − 1 = 8.

Homography and Action of a Projective Camera on a Crater Disk Quadric
If a crater rim is a planar conic, then its projection into an image is also a conic. The projective transformation from one plane to another is a homography, thus allowing the appearance of the projected crater ellipse to be computed analytically.
If C i is a 3D conic (see Eq. 38 and Table 1), we define its perspective projection into an image as A i = π(C i ), where π : P 3 \ {r M } → P 2 andr M is the camera location. Thus, we describe the conic locus for A i within the image as . We now introduce a coordinate system to make this projection explicit.
Consider a calibrated camera with selenographic position r M and attitude T M C that views lunar crater i. If ξ M is the selenographic position of a point on the rim of crater i, then its location in the camera frame is simply As discussed in Ref. [26], we can construct the pinhole camera by lettingx C ∝ ξ C , wherex C is the projected image plane coordinate of the 3D point ξ C . With a digital camera, however, we never observe a pointx C , but instead we observe the 2D pixel corresponding to that location in the image plane. Assuming the camera frame convention from Refs. [24] and [26], we place the +z axis out of the camera and along the optical axis, the +x axis to the image right and in the direction of increasing pixel column number, and the y axis completing the right-handed system. Furthermore, let the origin of the pixel u-v system be in the upper left-hand corner of the image, the u axis be to the right (increasing pixel column count), the v axis be down (increasing pixel row count), and integer values of [u, v] occurring at the pixel centers. With these conventions, the conversion between image plane coordinates (x C ) and digital image pixel coordinates (ū) is described by a simple affine transformation, or, more compactly,ū where the camera calibration matrix K is a function of the five calibration parameters: d x (and d y ) is the ratio of the focal length to pixel pitch in the x (or y) direction, α is the detector skewness, and [u p , v p ] is the coordinates of the principal point (where the optical axis intersects the image). We generally have excellent knowledge of K for spacecraft OPNAV applications from in-flight calibration with star field images [13,24]. Therefore, the 3D pointξ M projects to pixel coordinateū in the image according toū where P M C is the camera projection matrix for a specified absolute pose, The action of a projective camera on a quadric envelope is known to follow [51] where Q * i is from Eq. 40 and where A i (which may be computed from A * i ) describes the conic in the image plane tracing the apparent outline of the quadric. This allows for the analytic projection of the crater's disk quadric into its apparent ellipse in the image plane.
The result of Eq. 47 may also be viewed as a homography. Returning to Eq. 45, substitute the result from Eq. 33 forξ M , where we define H C i to be the 3 × 3 matrix describing the homography between the crater plane P i and the camera's image plane Since H C i is a homography, the conic locus and conic envelope may be analytically transformed to the image plane. Observing that we arrive at the result It follows, therefore, that The result of Eq. 56 is clearly the same as Eq. 47. We often find it convenient to make explicit the relative scales of A i and C i , which may be done by introducing the scale s i . For example,

Existence of Invariants from Projections of Crater Patterns
The central premise of this work is that invariant theory is the proper framework for describing, indexing, and matching patterns of lunar crater rims as seen in digital images. Recognizing that invariant theory will be unfamiliar to many space scientists and engineers, we find it worthwhile to develop the concept fully.
Much of the existing literature pursues an ad hoc approach for the construction of crater pattern descriptors. The lack of a formal framework has led to other authors proposing descriptors of varying dimension (ranging from two for a pair of craters [18,20] to 30 for a triad of craters [110]). However, for a given type of crater pattern there is a specified number of independent projective invariants. Descriptors with fewer elements than this do not fully exploit the perceptually salient information in the image. Descriptors with more elements than this increase complexity of the matching process without adding any perceptually salient information (the redundant invariants are all algebraic functions of the independent invariants). Thus, ad hoc schemes for developing crater pattern descriptors generally provide suboptimal performance, either in terms of descriptive power or computational complexity.

Preliminaries on Invariant Theory
In invariant theory we study functions on a space that remain unchanged under a group of symmetries of that space. Invariant theory originated in the 19th century with pioneering work by Cayley, Clebsch, Gordan, Sylvester, Hilbert, and others. Groups of particular interest were the general linear group GL n of invertible n × n matrices, the special linear group SL n of n × n matrices with determinant 1, orthogonal groups, and finite groups (see [35,36]). The space of binary forms of degree d (homogeneous polynomials in 2 variables of degree d) with the SL 2 -action described below was extensively studied in the 19th century [106]. If a(x, y) = a 0 x d + a 1 x d−1 y + · · · + a d y d is a binary form of degree d, then a matrix p q r s acts on it by where a 0 , . . . , a d are polynomial expressions in a 0 , a 1 , . . . , a d , p, q, r, s. For d = 2, the discriminant a 2 1 − 4a 0 a 2 is invariant under the SL 2 -action and every other polynomial invariant is a polynomial expression in the discriminant. Hilbert's papers [56,57] on invariant theory form the foundation of modern algebraic geometry in which one studies algebraic varieties. An algebraic variety is a set that is defined by polynomial equations, and an algebraic group is an algebraic variety that also has a group structure. Over the 20th century, invariant theory has been generalized to an action of arbitrary algebraic group G on an arbitrary algebraic variety V .

Polynomial Versus Rational Invariants
Traditionally the focus has been on polynomial invariants; i.e., polynomial functions that remain unchanged under the group symmetries. In practice, we generally desire a set of fundamental invariants. That is, we seek a set of invariants f 1 , f 2 , . . . , f r such that every other polynomial invariant g is a polynomial expression in the fundamental invariants: g = G(f 1 , f 2 , . . . , f r ) for some polynomial G(x 1 , . . . , x r ) in r variables. In the language of commutative algebra, the set of all polynomial functions on the variety V form a commutative ring, which we write as  [56,57] states that there is a finite list of fundamental invariants for groups such as GL n , SL n , and orthogonal groups. More generally, there is a finite system of fundamental invariants when the symmetry group is reductive [48,99]. For a precise definition of reductive, see Ref. [61]. Not all groups are reductive though, and typically, translation groups may not be reductive. There are examples where there is no finite set of fundamental invariants [98]. In many situations the smallest possible number of fundamental polynomial invariants is finite but extremely large compared to the dimension of the space.
In practical applications, such as computer vision and spacecraft optical navigation, it may not be feasible to find a fundamental system of invariants. One remedy is to work with separating invariants (see [34, §2.4]), and it is even simpler to work with rational invariants. Thus, instead of dealing with polynomial functions, we will consider rational functions.
Let C(V ) be the set of all rational functions on V , and C(V ) G be the set of all invariant rational functions on V . We assume that V is an irreducible variety so that C(V ) has the algebraic structure of a field-specifically, it is the quotient field of the ring C[V ]. We now see that C(V ) G is a subfield of C(V ). The reader should be warned that the quotient field of C[V ] G can be strictly smaller then C(V ) G . For example, it is possible that there are only constant polynomial invariants while there are non-constant rational invariants. This is another reason why it may be preferable to work with rational invariants instead of polynomial invariants.
Rational functions are not defined for all points in the space, but are defined for almost all points-which is good enough for practical crater identification. Invariant rational functions f 1 , f 2 , . . . , f r form a system of fundamental rational invariants if every other rational invariant g is a rational expression in f 1 , f 2 , . . . , f r ; i.e., of the form . , x r ) and G 2 (x 1 , . . . , x r ) are polynomials with coefficients in C and G 2 (f 1 , . . . , f r ) = 0. There always exists a finite system of fundamental rational invariants f 1 , f 2 , . . . , f r . In fact, it is possible to choose r ≤ 1 + dim V (see Remark 1).

Real Versus Complex Rational Invariants
Though the problem of spacecraft optical navigation (and computer vision, more generally) deals with geometric configurations over the real numbers, we will often work over the complex numbers. There is not a big difference. In the typical setup, we have a real algebraic variety V R that parametrizes certain geometric configurations. For example, V R could be a subvariety of the real projective space P m R defined by homogeneous polynomial equations with real coefficients. The same equations define also a projective variety V C ⊂ P m C in complex projective space if we view those equations over the complex numbers. We now assume that V R is Zariski dense in V C , which means that any complex polynomial that vanishes on V R also vanishes on V C . In this case, we can view the field R(V R ) of real rational functions on V R as a subfield of C(V C ). Moreover, the real and complex part of a rational function in C(V C ) are rational functions on V R . This gives a decomposition We also assume that G C is a complex algebraic group acting on the variety V C , such that the subgroup G R of real group elements acts on V R and G R is Zariski dense in G C . Then, because G R ⊂ G C is Zariski dense, a complex rational function is invariant under G C if and only if it is invariant under the group G R . Thus, we find that . . , f r )), then they also generate . . , f r , then the real and complex parts of f 1 , f 2 , . . . , f r generate the field R(V R ) G R of real rational invariants.

Counting Independent Rational Invariants
As we have seen, the proper counting of independent invariants is a recurring problem in the crater identification literature. Moreover, the usual method of counting invariants used in classical computer vision studies (e.g., a degrees-of-freedom approach) is not extensible to the problem at hand. We find it necessary, therefore, to develop a more formal framework. This framework is shown to reproduce known results for simple cases, before being extended to the problem of non-coplanar conics.
Suppose that L is a field and K is a subfield. If f 1 , f 2 , . . . , f r ∈ L, then K(f 1 , f 2 , . . . , f r ) is the set of all rational expressions in f 1 , . . . , f r with coefficients in K. This is also called the field that is generated by In spacecraft navigation, computer vision, and other applications it is often desirable to know the maximum number of functionally independent invariants. This relates to the notion of transcendence degree of a field extension. If f 1 , f 2 , . . . , f r lie in a field, then we say that f 1 , f 2 , . . . , f r are algebraically dependent over a subfield K if there exists a nonzero polynomial P with coefficients in K with P (f 1 , f 2 , . . . , f r ) = 0. The transcendence degree of L over the subfield K is the supremum over all r for which there exists f 1 , f 2 , . . . , f r ∈ L that are algebraically independent over K. Let trdeg(L/K) be the transcendence degree of L over K.
If L is also a subfield of another field M, then we have the following relation (see [ This shows that C(V )/K is algebraic; i.e., every element of C(V ) is algebraic over by the Theorem of the Primitive Element (see [72,Theorem 4.6]).
We are particularly interested in trdeg(C(V ) G /C) which is the maximal number of rational invariants that are algebraically independent (over C).
is the largest dimension of an orbit, then almost all orbits in V have dimension s and we call s the generic dimension of an orbit in V . The following results follows from Rosenlicht's Theorem [122].

Theorem 1 If s is the dimension of a generic orbit in V , then the maximum number of algebraically independent invariants is equal to
Proof Rosenlicht proved in [122] that there exist a G -stable nonempty (Zariski) open subset U ⊆ V that has a geometric quotient. A geometric quotient is an algebraic variety U /G together with a surjective morphism π : U → U /G such that the fibers of π are exactly all G -orbits and π also has additional properties such as

Rational Invariants for Conics In P 2
The result of Theorem 1 is straightforwardly applied to the problem of counting the rational invariants for coplanar conics (i.e. a d-tuple of conics in P 2 ). We know from past results that two coplanar conics have two invariants [41,96,115], that three coplanar conics have seven invariants [116], and that a d-tuple of conics have 5d − 8 invariants [53]. This is now briefly shown, before moving on to the more nuanced problem of non-coplanar conics.
Therefore, as an illustration, we can count the number of independent rational invariants for d conics in P 2 . A conic in P 2 corresponds to a nonzero quadratic homogeneous polynomial in 3 variables, up to scalar (see Eq. 9). So the variety of all conics can be identified with P 5 because such a polynomial has 6 = 5 + 1 coefficients. Let V be the variety of d-tuples of conics in P 2 , Then V ∼ = (P 5 ) d has dimension 5d. For G we take the group PGL 3 which has dimension 8. For d = 1 the group PGL 3 acts transitively on all nondegenerate quadratic forms. This means that PGL 3 has a Zariski dense orbit in V = P 5 , so the dimension of a generic orbit is s = 5 and the number of independent rational invariants is dim V − s = 5 − 5 = 0. For d = 2, we can explicitly compute the stabilizer of a pair of conics (e.g., x 2 + 2y 2 − z 2 = 0 and 2x 2 + y 2 − z 2 = 0) and note that it is finite. This implies that the dimension of a generic orbit is equal to dim G = 8. This, in turn, implies that the dimension of a generic orbit is equal to 8 for all d ≥ 2 and the number of algebraically independent invariants is dim V − s = 5d − 8.

Rational Invariants for Conics in P 3
A method for the robust identification of non-coplanar crater patterns has eluded spacecraft navigators since autonomous crater-based navigation was first explored over 25 years ago. The result is that most existing algorithms constrain the problem to local (nearly coplanar) crater patterns or to only nadir pointing images. To solve this problem requires us to first understand conics in P 3 .

Counting independent rational invariants for conics in P 3
Let Z be the 8-dimensional variety introduced at the end of Section 4.1.2 parameterizing conics in P 3 , and let V = Z d be the space of d conics in P 3 . The 15-dimensional group G = PGL 4 acts on P 3 , Z , and V . If d = 1 then the group acts transitively on all nondegenerate conics so that V = Z has a dense orbit and there are no non-constant rational invariants. Thus, we search for invariants when d ≥ 2.
Suppose that d = 2 and consider the pair of conics C 1 and C 2 . If the conics are generic enough, then the stabilizer of the pair (C 1 , C 2 ) is finite. To see this, suppose that g ∈ PGL 4 lies in the connected component of the identity in the stabilizer of (C 1 , C 2 ). Then g must also fix the planes P 1 and P 2 through C 1 and C 2 respectively. Further, g must also fix the intersections C 2 ∩ P 1 = {c (1) 2 ,c (2) 2 } and C 1 ∩ P 2 = {c (1) 1 ,c (2) 1 } and because g lies in the connected component of the identity, it fixes all the pointsc 2 on the line P 1 ∩ P 2 individually. There are two points r (1) 1 ,r (2) 1 on C 1 such that the tangent lines to C 1 atr 1 in the plane P 1 there are no three points on the line (Fig. 9), because they all lie on the same conic C 1 . This implies that g must fix all the points in the plane P 1 . Similarly it must fix all the points in the plane P 2 . From this follows that g is the identity element in PGL 4 . Thus, the dimension of a general orbit is equal to dim PGL 4 = 15, and the same is true for all d ≥ 2. It follows, therefore, that the maximal number of algebraically independent rational invariants for d ≥ 2 conics is dim This result agrees with the past literature that has identified one independent rational invariant for a pair of conics in P 3 [65,114]. Our result, however, is more general and prepares us to study which invariants (if any) may be recovered from the projection of conics into an image.

Calculating Invariants from Projections
In Section 5.3.1, we found the invariants for a d-tuple of conics in P 3 . These, however, are not the invariants we need because we do not measure the crater rims directly in P 3 . Instead, we image the lunar surface with a camera. Thus, we are in search of rational invariants that may be computed from only the projected crater rims that we see in an image. Such invariants do not always exist.

Invariants from Projections and the Action of PGL 4
We begin by discussing the general setup for computing invariants from projections. To do this, we consider models M in C 3 space that are parameterized by some variety 1 ,c (1) 1 ,c (2) 1 all lie on C 1 no three of these points may be on the same line V . For example, V can consist of all d-tuples of points in 3D, or all d-tuples of conics. We also have camera positionō ∈ C 3 and a perspective projection π : C 3 \ {o} → P 2 . We are interested in rational functions on V (functions depending on the model M) that can be computed from the projected image π(M) in a way that is independent on the camera positionō. We can extend C 3 to the projective space P 3 and also extend the camera projection to a map P 3 \ {ō} → P 2 , whereō is equal to o ∈ C 3 ⊂ P 3 , but viewed in P 3 . The group PGL 4 acts on P 3 and contains the group of affine transformations of C 3 . The group PGL 4 does not act on C 3 , but it acts on the field of rational functions on C 3 , because C 3 and P 3 have the same rational functions. We assume that the class of models parameterized by V is closed under the action of PGL 4 (for example, PGL 4 takes conics to conics). Then PGL 4 also acts on the rational functions on V .
Suppose that f is a rational function on V such that f (M) can be computed from the projected image π(M) in such a way that it does not depend on the choice of the positionō of the camera. So there exists a function h such that f (M) = h(π(M)) for every choice ofō. Let Gō be the image of the group in PGL 4 . An element g ∈ Gō fixesō and all the lines throughō, so we have So f is invariant under the action of Gō. Since this is true for every choice ofō, we see that f is invariant under the group G generated by all Gō,ō ∈ P 3 . One can verify that G = PGL 4 . To see this, for example, we observe that the group G is closed under conjugation with elements in PGL 4 . This implies that G is a non-trivial normal subgroup of PGL 4 . Since PGL 4 is known to be a simple group, we get G = PGL 4 . This shows that any rational function that can be computed from a camera projection in a way that is independent on the choice of the position of the camera must be invariant under the action of PGL 4 .

Invariants from Projections of Points in P 3
The simplest example of invariants from projections is the case of a d-tuple of points in P 3 . It is well known that d arbitrarily placed 3D points possess no invariants that one may compute from their projection in an image [16,28]. Thus, a cloud of random 3D points cannot be indexed for recognition by a pose invariant descriptor, which is a critically important fact. We reproduce this known result here in our more formal framework. In doing so, we develop some of the ideas and tools necessary for the case of conics in P 3 but within the context of a simpler (and familiar) example. Therefore, for example, suppose we have d points in P 3 . For d ≥ 5 there are 3d −15 independent invariants for the action of PGL 4 (because for d ≥ 5, the generic stabilizer in (P 3 ) d is finite). Perspective projection for some fixed cameragives d points in P 2 . Unfortunately, however, no non-constant rational invariant can be computed from just the image. This may seem counter-intuitive, and we now show why this is the case.
Suppose we fix a model M ∈ (P 3 ) d consisting of d points. The image may change in appearance with varying camera viewpoint, with the projection being governed by the action of PGL 4 . The variety of the possible images from the model as we move the camera cuts out a variety U ⊆ (P 2 ) d of dimension at most dim PGL 4 = 15. For d ≥ 8, the dimension of U is strictly smaller than the dimension of (P 2 ) d and one might expect to have at least 2d − 15 independent invariants that can be computed from the images, but this is wrong as we will see. The problem is that Theorem 1 does not apply. The variety U of possible images is not an orbit for any group action on (P 2 ) d : the group PGL 4 does not act on (P 2 ) d , and U may be bigger than any PGL 3 orbit.
To proceed, suppose that V = (P 3 ) d and Y = (P 2 ) d ,ō ∈ P 3 is the position of the camera, and π : V → Y is the projection. (The map π is only defined for elements in (P 3 We are interested in rational invariants that can be computed from the image in Y . Such invariants are exactly elements in the intersection field C(Y ) ∩ C(V ) G . We get the following diagram of field extensions with their transcendence degrees: We will see that C(Y ) ∩ C(V ) G = C, which means that that there are no nonconstant rational invariants that can be computed from the image.
To understand what is happening, we introduce an equivalence relation ∼ on (P 2 ) d , the variety of possible images of d points in P 2 . We say I 1 ∼ I 2 is true when there exists a model M such that I 1 , I 2 both appear as images of M (but possible from different camera positions). The relation ∼ is not an equality relation. It satisfies the reflexivity axiom (I ∼ I) and the symmetry axiom (I 1 ∼ I 2 ⇔ I 2 ∼ I 1 ). However, it does not satisfy the transitivity axiom (if I 1 ∼ I 2 and I 2 ∼ I 3 , then I 1 ∼ I 3 ). If π and π are the projections with respect to different camera positions/orientations, then π is obtained from π by a projective linear transformation g, so that π (M) = π(g · M) for any model M. If I ∼ I , then we have I = π(M) and I = π (M) for some model and some camera projections π, π . This means that h(I) = h(π(M)) = f (M) = f (g · M) = f (π(g · M)) = h(I ).
Let ≡ be the equivalence relation generated by ∼. So we say that I ≡ I if and only if there is a finite sequence of images I = I 0 , I 1 , I 2 , . . . , I r = I such that I 0 ∼ I 1 , I 1 ∼ I 2 , . . . , I r−1 ∼ I r . If I ≡ I and I 0 , I 1 , . . . , I r are as above, then f (I) = f (I 0 ) = f (I 1 ) = · · · = f (I r ) = f (I ).

Invariants from Projections of Conics in P 3
The corollary to a d-tuple of arbitrary points in P 3 is a d-tuple of arbitrarily placed conics in P 3 . We will show that there are no non-constant invariants for d conics in P 3 that can be computed from a perspective projection. The argument is similar to the argument in Section 5.4.2 using an equivalence relation on d-tuples of conics in P 2 . We use an explicit geometric construction to show that all d-tuples are equivalent.
As before, let Z be the variety of conics in P 3 defined at the end of Section 4.1.2. Then, let V = Z d be the variety of d-tuples of conics in P 3 and Y = (P 5 ) d be the variety of d-tuples of conics in P 2 . We define a relation ∼ on Y by (A 1 , A 2 , . . . , A d ) ∼ (A 1 , A 2 , . . . , A d ) if there exists two camera projections π and π and conics (C 1 , Let ≡ be the equivalence relation generated by ∼. Suppose that B 1 , A 1 , A 2 , . . . , A d ⊆ P 2 are conics. Let (1) , (2) be two lines in P 2 that are tangent to both A 1 and B 1 with A 1 , B 1 within the same region of P 2 \ ( (1) ∪ (2) ) (see Fig. 11).
We denote the point at which (j ) is tangent to A 1 (respectively B 1 ) byā . Also, letq be the intersection point of (1) and (2) . Choose a pointō ∈ P 3 (the center of the camera) and let π : P 3 \ {ō} → P 2 be the camera projection. Let X i = π −1 (A i ) ∪ {ō} be the cone with topō that corresponds to A i ⊂ P 2 . Also, define Y 1 = π −1 (B 1 ) ∪ {ō}. The lines (j ) correspond to the plane L (j ) = π −1 ( (j ) ) ∪ {ō} for j = 1, 2. The planes L (1) and L (2) intersect in the line Q = π −1 (q) ∪ {ō}, as shown in Fig. 12. We choose another pointō on the line Q and π : P 3 \ {ō } → P 2 is the camera perspective projection with respect to the camera positionō . We choose some random plane P in P 3 . Then Fig. 11 Visualization of geometric quantities used to show the non-existence of invariants from the projection of arbitrary conics in P 3 . The conics and lines in this graphic all lie within the image plane. The backprojection of these planar quantities into P 3 form quadric cones (for A 1 and B 1 ) or planes (for (1) and (2) )

Fig. 12
Visualization of geometric quantities used to show the non-existence of invariants from the projection of arbitrary conics in P 3 . The lines T (1) and S (1) lie in plane L (1) (lines are white, plane is gray). The lines T (2) and S (2) lie in plane L (2) C j = P ∩ X j is a conic in P 3 and π(C j ) = A j . Let X j be the cone with topō through C j . By construction, we have π(C 1 , We will construct a conic D 1 in P 3 with π(D 1 ) = B 1 and π (D 1 ) = A 1 . Note that X 1 is a cone with topō that is tangent to the planes L (1) and L (2) . The intersection 1 )∪{ō} is a line throughō, and T (j ) = X 1 ∩L (j ) is a line throughō for j = 1, 2. Letr (j ) be the intersection point of the lines S (j ) and T (j ) in the plane L (j ) for j = 1, 2. Then we have π(r (j ) ) =b (j ) 1 . Moreover we choose a pointr (3) =r (1) ,r (2) that lies on Y 1 ∩ X 1 . Let R be the plane throughr (1) ,r (2) ,r (3) . Define D 1 = R ∩ X 1 . Then D 1 is tangent to the planes L (1) and L (2) at the points r (1) andr (2) respectively. Now π(D 1 ) and B 1 both are tangent to (1) and (2) at the pointsb (1) = π(r (1) ) andb (2) = π(r (2) ) respectively, and both contain the point π(r (3) ). It follows that π(D 1 ) = B 1 because a conic in P 2 is determined by 3 points and the tangent lines at 2 of the points. One would expect the intersection Y 1 ∩ X 1 of two quadratic surfaces to be a curve of degree 4. In this case, the curve is reducible, and a union of two conics, and D 1 is one of them.
Repeating the argument shows that (A 1 , A 2 , . . . , A d ) ≡ (B 1 , B 2 , . . . , B d ) for all conics A 1 , A 2 , . . . , A d , B 1 , B 2 , . . . , B d . This implies that there is no nonconstant rational invariant for d conics in P 3 that can be computed from a projection.
That no rational invariants exist for the projection of d arbitrary conics is an important finding (which we believe to be a novel result). This clearly precludes the use of invariants for lost-in-space crater identification about an arbitrarily shaped body. Fortunately, celestial bodies large enough to exhibit substantial cratering do not have an arbitrary shape. The surface of these bodies can usually be modeled regionally-if not globally-as a nondegenerate quadric surface. Specifically, we note that planets, moons, and dwarf planets are generally ellipsoidal in global shape [90] (this is certainly the case for the Moon), while the smaller minor planets are often globally irregular but regionally ellipsoidal. Comets nuclei often show no such structure. Further, many large bodies (e.g., the Moon) appear nearly planar over sufficiently small portions of the surface. We know from before that invariants exist for coplanar conics (Section 5.2). We now show that invariants also exist for non-coplanar conics lying on a quadratic surface (a quadric).

Counting Invariants for Conics on a Nondegenerate Quadric Surface
As we have seen, there are no nontrivial invariants for d-tuples of arbitrary conics in 3D space that can be computed from a camera image. In the case of craters on the Moon, the conics we are interested in are not arbitrary, but lie on the surface of the Moon. Thus, instead of considering d-tuples of arbitrary conics, we consider the variety V of d-tuples of conics for which there exists a nondegenerate quadric surface that contains all of them. We find projective invariants to exist for d ≥ 3 conics lying on the same nondegenerate quadric surface.

Pairs of Conics on a Nondegenerate Quadric Surface
First consider the case d = 2. Suppose that C 1 and C 2 are two conics that lie in the planes given by the linear equations f 1 = 0 and f 2 = 0 respectively. Then C 1 and C 2 lie in the degenerate quadric surface f 1 f 2 = 0. Suppose both conics also lie on a nondegenerate surface defined by a quadratic equation g = 0. This surface is not unique, because for every t, the surface defined by g + tf 1 f 2 = 0 contains both conics. If two conics C 1 and C 2 lie on a nondegenerate surface defined by g = 0, then C 1 is defined by f 1 = g = 0 and C 2 is defined by f 2 = g = 0. The intersection C 1 ∩ C 2 is defined by f 1 = f 2 = g = 0 and consists of two (possibly complex) points. On the other hand, if C 1 ∩ C 2 consists of two points, one can show that both conics lie on a nondegenerate quadric surface. To find a pair (C 1 , C 2 ) in V , we can choose the conic C 1 arbitrarily in the 8-dimensional variety of conics, we choose the plane P 2 that contains C 2 in the 3-dimensional space of hyperplanes, Now the conic C 2 is an element of the 5-dimensional variety of conics in P 2 , but there are two constraints because the conic has to go through the two points of C 1 ∩ P 2 . So the dimension of the variety V is 8 + 3 + (5 − 2) = 14. The 15-dimensional group PGL 4 acts on V . Let H be the connected component of the identity in the stabilizer of (C 1 , C 2 ) ∈ V . The conics C 1 and C 2 intersect in two (possibly complex) pointsc 1 andc 2 . So g will fix the pointsc 1 andc 2 . If L is the line throughc 1 andc 2 then g will map L to itself. Ifc 3 is a 3rd point on L then g may mapc 3 to another point of L. However, if g also fixes the pointc 3 then a similar argument as in Section 5.3.1 shows that g must fix the planes P 1 and P 2 pointwise, and must be the identity. This shows that the stabilizer of the pair (C 1 , C 2 ) is at most 1-dimensional. Since the generic stabilizer has dimension at most 1, the dimension of a generic orbit is at least 15 − 1 = 14. Since dim V = 14, this implies that there are no rational invariants for a pair of conics lying on a nondegenerate quadric surface.

Many (d ≥ 3) Conics on a Nondegenerate Quadric Surface
Let us now assume that d ≥ 3. If a quadric surface S contains 3 distinct conics C 1 , C 2 , C 3 then S is the unique quadric surface that contains these three conics. To see this, suppose that S is another quadric surface through the 3 conics. Suppose that S and S are defined by the equations h = 0 and h = 0 respectively, where h and h are homogeneous quadratic polynomials in 4 variables. Let P i be the plane through C i given by f i = 0 where f i is a linear function. Letā ∈ P 3 . We can multiply h and h with nonzero scalars such that (h − h )(ā) = 0. In the plane P i , the restriction of the To parametrize d conics that lie on a quadric surface, we can first choose the surface that is given by a quadratic polynomial in 4 variables up to a scalar. Such a polynomial has 10 coefficients, so the quadric surface is determined by 10 − 1 = 9 parameters. Now, each of the conics is determined by a hyperplane section of the quadric surface. Hyperplanes are parameterized by points in P 3 . So the variety V of d-tuples of conics that lie on a common quadric surface has dimension 9 + 3d. The stabilizer of a generic point in V is finite, so the dimension of a generic orbit is 15 = dim PGL 4 . So the number of independent rational invariants is (9+3d)−15 = 3d −6 for d ≥ 3.

Computing Invariants from Crater Rims in an Image
The results of Section 5 established the existence of invariants for conics lying on either a plane or a nondegenerate quadric surface. While the coordinate-free approach used in the previous section is a powerful tool for studying such invariants, we must ultimately impose a coordinate system for practical computation of these invariants from data. Doing so is straightforward, but does require some additional mathematical machinery. The details are now discussed.

Pair of Arbitrary Non-Coplanar Conics
Two arbitrary 3D conics do not generally intersect one another, even over the complex numbers. This may be seen through the intersection of surfaces (illustrated in Fig. 13). The line formed by the intersection of the two conic planes intersects each conic at two places, producing four intersection points. The cross-ratio of these Fig. 13 Two arbitrarily placed non-coplanar conics do not intersect one another. The common 3D line L ij formed by the intersection of their planes intersects each of the two conics in two places, creating four (usually distinct) points. The two intersection points for a particular conic are complex valued if L ij does not physically intersect that conic. The four colinear intersection points (possibly complex valued) may be used to form a cross-ratio, which is the single unique 3D invariant for a pair of 3D conics four colinear points is the unique 3D invariant for a pair of non-coplanar conics. This was observed in [87] and [114], which we now rederive by other means as we build towards a methods for computing projective invariants. Further, the results of Section 5.3.1 tell us that this is the only 3D invariant for a pair of non-coplanar conics. We also show that this 3D invariant does not lead to a useful projective invariant that may be constructed from an image of these two conics.
Suppose we have crater i described by the 3D conic C i that lies in plane P i . Let the Moon-centered quadric cone of crater i be given by X i as described in Eq. 37. By construction, this cone must pass through the conic locus C i , and we interpret C i as the conic section from Eq. 38. Now, suppose we have two craters: crater i and crater j . We compute their intersection by substitution of Eq. 38 as, Define the 3D line L ij as the intersection of planes P i and P j , such that The intersection of a line with a quadric cone, e.g. L ij ∩ X i , produces two points. Define the intersection points of this line with crater i as {c (1) i , c (2) i }. It is easy to see, as illustrated in Fig. 13, that these intersection points are not generally the same for two arbitrary conics. Therefore, Observe that the four points {c (1) i , c (2) i , c (1) j , c (2) j } are colinear since they all lie on the line L ij . They are also distinct, meaning their cross ratio ρ ij is a 3D invariant of the conic pair [114] ρ ij = Cr(c (1) i , c (2) i , c (1) j , c (2) j ) The difficulty with this 3D invariant is that it is not recoverable from a single image of the two conics. The two 3D conics will project into two coplanar conics in the image plane. The image conics have four (often complex valued) intersection points, which are not necessarily colinear. In the absence of other constraints, these intersection points in the image plane are not related to the 3D line L ij or its intersection points {c (1) i , c (2) i , c (1) j , c (2) j } with the two 3D conics. This is because the two conics do not actually intersect one another. That no projective invariant exists from a pair of non-coplanar conics agrees with the findings of Section 5.4.3. The non-existence of projective invariants for three or more arbitrarily placed non-coplanar conics follows by similar arguments.
In the special case where the conics do intersect one another, their intersection is preserved when imaged from an arbitrary pose by a projective camera. One of the most flexible ways to ensure two 3D conics intersect over the complex numbers is to constrain them to lie on a nondegenerate quadric surface (e.g., sphere, ellipsoid). This is now shown.

Pair of Non-Coplanar Conics on a Nondegenerate Quadric Surface
Consider lunar craters modeled as conics lying on a nondegenerate quadric surface S. Under this assumption, we may model the surface of Moon by the quadric locus parameterized by Q M , We observe here that the quadric surface need not be a global surface fit, but only a good approximation of the portion of the lunar surface where the craters reside. Indeed, some crater clusters may use one quadric surface while other crater clusters use a different quadric surface. This choice in no way affects the theoretical developments that follow. Now, consider crater i lying on plane P i . For the 3D rim of a crater i to also be on the surface of the Moon, it must lie on the intersection of this plane with the S. Therefore, the 3D conic locus for crater i must be the intersection which produces a circular crater if S is a sphere. An ellipsoidal (non-spherical) body would generally produce an elliptical crater. The discussion that follows holds for both spherical and ellipsoidal bodies. Further, as in Section 6.1.1 and described in Eq. 38, the crater may also be viewed as the conic section produced by intersecting plane P i with the Moon-centered quadric cone X i . However, since the crater must lie on both X i and S, it is clear that we can also write the 3D crater conic as, Therefore, we see that Suppose we have two craters: C i and C j . Because C i and C j are formed by the intersections of planes P i and P j with the quadric surface S, we observe that, where L ij is the 3D line formed by the intersection of the two crater planes (see Eq. 60). This is illustrated in Fig. 14.
If the two craters physically intersect, the line L ij pierces S and the intersection points are real. If the craters do not physically intersect, the intersection occurs over the complex numbers. Lets (1) ij ands (2) ij be the two intersection points {s (1) ij ,s (2) Therefore, define P ij as the plane passing through the center of the Moon and the line L ij (or, equivalently passing through the center of the Moon,s ij , ands ij ). Clearly, the numerical stability of computing P ij in this manner is poor if the craters are close to one another and P i and P j are nearly coplanar. Better numeric stability may be achieved by taking a different approach.
The conic intersection from Eq. 68 may also be written in terms of the quadric cones instead of the planes, Without constraints, the intersection of two cones X i ∩X j is a set of four lines passing through the center of the Moon that pierce S at eight different points. Four of the points will be on the wrong side and come from conics that do not actually exist; these can be ignored. When the conics lie on a nondegenerate quadric surface, we 14 Two non-coplanar conics C i and C j lying on a nondegenerate quadric surface must intersect in two points (which happen to be complex valued in this example). The two intersection points lie on the line L ij formed by the intersection of the two crater planes. The intersection points also lie in the plane P ij formed by the join of L ij and the body center find the other four intersection points are a set of repeated intersections. To see this, combine the above results to find, ij ,s (2) ij } The advantage of using the quadric cone intersections to describes (1) ij ands (2) ij is that this remains well defined and more numerically stable as the craters become close and their planes become nearly parallel. This is especially important when viewing craters at a local level, where craters are nearly coplanar.
Since the pointss (1) ij ands (2) ij lie on C i , the projection of these points lie on the projection of C i (which we called A i ). The same holds for C j , and its projection A j . Therefore,the projection ofs (1) ij ands (2) ij is recoverable as two of the intersection points of the image conics. That is, where P M C is the projection matrix from Eq. 46 and the pixel coordinateū ij is the projection ofs (k) ij into the image. The challenge, therefore, is to determine which of the four intersection points of A i and A j are the projection ofs (1) ij ands (2) ij . There are six possible combinations of these four points, corresponding to six possibilities for the projection of line L ij into the line ij . We now discuss how to uniquely compute ij .

Consider the two conics A i and A j described by the conic locus matrices A i and A j . We recall thatū
It follows that we can form a pencil of conics parameterized by the scalar ratio λ : μ that also pass through the same four intersection points as A i and A i , If we force the matrix λA i + μA j to be a degenerate conic, then it becomes two lines (or a double line) and we may find its intersection with the two conics. We find six possibilities, one of which is the line ij that we seek. Letting μ = 1 for easy computation, we force the conic to be degenerate by setting the determinant to zero, Since both A i and A j are full rank and well-conditioned, we may rewrite this as which is a simple 3 × 3 eigenvalue problem. When there is no actual intersection of the two conics in the image, only one of the resulting degenerate conics will produce real valued lines [117]. The specific line we seek is ij , which is the projection of L ij . As the intersection of two non-coplanar planes, L ij is real valued by construction. Thus, its projection ij is also real valued. Therefore we always seek the eigenvalue leading to a degenerate conic of real valued lines, discarding the four complex valued line possibilities. It will soon become apparent which eigenvalue to choose. Therefore, suppose we choose one of the eigenvalues to construct the degenerate conic formed by where we observe that B ij has rank 2. The remaining steps are now to find the lines described by B ij , determine if they are real valued, and then determine which one is ij . A degenerate conic made of two lines, g and h, is defined by the symmetric matrix where our task is now to find g and h given B ij . There are a variety of ways to accomplish this task. Our specific approach is similar to the framework outlined in [117]. Therefore, proceeding in this way, we first recall that the intersection of the two lines in homogeneous coordinates is computed as [51] Further, observing that [z×] = gh T − hg T (80) we find that that is a rank one matrix with columns proportional to g and rows proportional to h. We must now find z.
To do this, a quick calculation will confirm that Or, in terms of the original ellipses, The diagonals of −B * ij are the squares of the elements ofz. Thus, defining the columns and elements of B * ij as any column of B * ij may be scaled to findz according tō where any k ∈ {1, 2, 3} with b kk = 0 will work. Best numerical performance is achieved by selecting the value for k that yields max k b * kk . If the ellipses do not intersect, only one eigenvalue from Eq. 76 will produce a real valued estimate of z. The eigenvalue to choose is the one that makes the diagonal of B * ij negative. Therefore, withz known, compute D we can then compute g from any non-zero column of D and h from any non-zero row of D. In general, it is best to find the element of D with the largest absolute value and pick the corresponding row and column for g and h.
Given the two lines g and h, one of these corresponds to ij . The other does not. Each of the lines g and h divide P 2 into two regions. Since the line we seek comes from the projection of L ij = P i ∩ P j , the two conics must not lie in the same region. That is, we seek a line that passes between the two conics. Only one line will satisfy this condition and this is the line we choose for ij .
If we were to compute the intersection of ij with either A i or A j we would obtain two points that are the projections (1) ij ands (2) ij . Fortunately such a computation is not necessary since we work directly with ij in subsequent discussions.
Briefly, we observe that C i and the two pointss (1) ij ands (2) ij all lie in the plane P i . While a projective invariant exists for two points and a conic (all coplanar) [41], the points must not lie on the conic. Sinces (1) ij ands (2) ij lie on C i by construction, we seek an alternative way of constructing an invariant.
Consider instead the 2D line ij formed by the projection of 3D line L ij (or, equivalently, by the join of image pointsū (1) ij andū (2) ij ). While there is not an invariant for a single line and a conic, there is an invariant for two lines and a conic [123]. This motivates the study of invariants for a triad of craters.

Triad of Non-Coplanar Conics on a Nondegenerate Quadric Surface
Suppose a projective camera observes three craters: C i , C j , and C k . The intersection of the corresponding planes (P i , P j , and P k ) produces the three 3D lines L ij , L ik , and L jk . Three planes intersect in a point, which is also the location where the 3D line formed by two planes intersects the third plane. This point, defined as a ij k , is the apex of the pyramid formed by the three planes, Now, let the three 3D craters C i , C j , and C k project to 2D conics A i , A j , and A k in the image. For any given pair of 3D conics (e.g., C i and C j ), the 3D line through their intersection points (e.g., L ij ) is coplanar with both of these 3D conics and it is possible to recover its projection (e.g., ij ) from just the projected image conics (e.g., A i and A j ). Thus, for a triad of observed craters in an image, we may compute the three lines ij , ik , and jk . Since C i , L ij , L ik are all coplanar, an invariant exists in the image using A i , ij , and ik . The same is true for craters C j and C k .
We may show that such an invariant exists by use of a cross ratio, and this invariant may be computed directly by use of a Cayley-Klein metric. The usual cross-ratio applies to four points on a line. However, by the duality of points and lines in P 2 , we can also form a cross-ratio of four lines passing through a point. Therefore, as shown in Fig. 15, consider a conic A i . Let the two lines ij and ik be described by the 3 × 1 vectors ij and ik , and let these two lines intersect at the pointp = ij × ik . Using the pole-polar relation for a conic, we may find the two lines fromp that are tangent to the conic A i , which we call w 1 and w 2 . Since these four lines go through the point p, we may form a cross ratio that is a projective invariant We can now form a classical Cayley-Klein metric as which, from hyperbolic geometry, is known to be equivalent to where the adjugate matrix A * i describes the conic envelope of A i . To understand the need for the absolute value in the numerator, recall from Section 4 that ij and − ij describe the same line (same for ik and − ik ). Likewise, A * i and −A * i describe the same conic envelope. Thus, since the appropriate choice of sign is not always clear a priori, the absolute value in the numerator is a simple way to robustly remove the sign ambiguity and ensure that cosh d( ij , ik ) ≥ 1 for all d( ij , ik ). Absolute values are not necessary in the denominator since this is quadratic in ij , ik , and A * i . We may also write, therefore, Since the cross ratio from Eq. 88 is a projective invariant, it follows immediately that the Cayley-Klein metric from Eq. 89 and its hyperbolic cosine in Eqs. 90 and 91 are also projective invariants. At first, the Cayley-Klein development may seem of purely academic interest. However, it is of critical importance when building a searchable index for lunar crater identification. It is only through this relation that we fully understand why the metric for α 2 i from Eq. 91 (the typical invariant found in the literature for a coplanar set of a conic and two lines [123]) loses descriptiveness for small distances d( ij , ik ). The spacecraft navigator will immediately see that this is analogous to how small interstar angles θ lose descriptiveness when cataloged as cos θ instead of by θ (i.e., we need to tighten the matching tolerance as θ → 0 when using a cos θ index). Thus, we prefer to index on d instead of α 2 .
We usually choose to only consider craters that do not physically intersect (i.e. do not overlap), since many CDAs have difficulty providing high-quality crater rim fits for intersecting craters. Moreover, if restrict ourselves to only non-overlapping craters, then we are guaranteed that the cross ratio is positive. This implies that d ≥ 0 and is real, and that the denominator on the right-hand side of Eq. 90 is also real. Since d ≥ 0, we know that Therefore, we may relate the cross-ratio, the Cayley-Klein metric ( ij , ik ), and the invariant α i according to where we choose to construct the invariant to index with the Cayley-Klein metric, An identical procedure produces corresponding invariants for craters j and k We know from our review of the literature that many previous attempts at crater identification have accidentally employed feature descriptors with elements that are not independent. Thus, before proceeding further, it is necessary to ensure that the invariants J i , J j , and J k are independent. Since we know from Section 5.5 that a dtuple of conics on a quadric surface has 3d − 6 invariants (leading to 9-6=3 invariants for this case), it follows that if the three invariants J i , J j , and J k are independent then we have found all the algebraically independent invariants that exist. This is now shown.
We observe that the invariants α 2 i , α 2 j , α 2 k (and hence also α i , α j , α k and J i , J j , J k ) are algebraically independent. We calculate these invariants for a particular model, since it suffices to show that the invariants are independent when restricted to special subclasses of models. We assume that the lunar surface is the sphere x 2 +y 2 +z 2 = 1, and that the three craters are the circles defined by intersecting the sphere with the planes given by x = t 1 , y = t 2 and z = t 3 respectively. Within the plane z = t 3 we have the crater given by x 2 + y 2 = 1 − t 2 3 and the lines x = t 1 and y = t 2 . So we get Thus, we may compute α 2 and, by symmetry, we also find It is possible to verify that the Jacobi matrix is invertible for some choice of t 1 , t 2 , t 3 (e.g., t 1 = t 2 = t 3 = 1/2). This implies that α 2 1 , α 2 2 , α 2 3 (and, therefore, J i , J j , J k ) are algebraically independent [75]. Despite the rather long development up to this point, the practical implementation of this method for computing the invariants J i , J j , J k is straightforward and is summarized in Algorithm 1. The reader should note that the invariants J i , J j , J k are computed directly from the measured image conics A 1 , A 2 , A 3 without requiring any knowledge of the camera position or orientation. This permits lost-in-space crater identification.
The utility of this framework is now briefly demonstrated. Consider a regional crater pattern such as the one shown in Fig. 16, where there is notable curvature of the Moon. This example has about 18.1 deg between the surface normal of crater A j and A k , thus necessitating the consideration of non-coplanar invariants. Two substantially different views of the same crater pattern are shown. The image on the left is nearly nadir pointing and the image on the right is pointed over 30 deg off nadir. As can be seen, the non-coplanar invariants J 1 , J 2 , and J 3 are identical in both images since Fig. 16 Example of invariants J i , J j , and J k (blue) for the same triad of craters (A i , A j , and A k ) seen from two different vantage points. The invariants remain unchanged under perspective projection from any viewing geometry. The triangle interior angles (yellow) are not a projective invariant. Synthetic images of the lunar surface were produced using PANGU Planet Surface Simulation Software developed by the Space Technology Centre at the University of Dundee, Scotland [111] they are formal projective invariants. This makes evident the power of such invariants for crater pattern recognition.
It should be stressed that only specific items in Fig. 16 are coplanar, with the vast majority of items being non-coplanar. Consequently general combinations of lines or crater ellipses from left image cannot be transformed to the right image with a common homography.
Also shown in Fig. 16 are the crater triad interior angles (yellow). We highlight these here because crater triangle interior angles are often proposed as descriptors for lost-in-space crater pattern recognition, e.g. [49]. It is clear, however, that the interior angles are not projective invariants and do not provide a robust means of pattern indexing if off-nadir viewing is possible.

Geometry of Invariants for Coplanar Conics
When viewing only a very small portion of the lunar surface, the observed craters are often nearly coplanar. Smaller craters are also more elliptical. Thus, patterns of small craters close to one another are better described using a triad of coplanar craters. Here, we revisit the problem of invariants of two and three coplanar conics using the framework of Semple and Kneebone and [128] and of Quan [115,116]

Pairs of Coplanar Conics
A pair of coplanar conics has two projective invariants [41]. We saw in Section 6.1.1 that two arbitrarily placed 3D conics have no intersections and in Section 6.1.2 that two 3D conics on a nondegenerate quadric surface have two intersections. We observe now that two coplanar conics always have four intersections. These four intersection points (possibly over the complex numbers, possibly repeated) persist under the action of a projective camera, which permits easy computation of two invariants.
Four specified points on a conic may be used to construct a cross-ratio, which is a projective invariant. It may initially be unclear why this is the case since the four points around the ellipse are not collinear, thus a brief explanation is warranted. As discussed earlier, using the duality of points and lines in P 2 , we can also form a projective invariant by the cross-ratio of four lines passing through a point. Chasles' theorem [128] states that the cross-ratio is a constant for the pencil of four lines from four points on a conic to a fifth point also on the conic (Fig. 17). Therefore, we may directly form two invariants: the first by finding the cross ratio of the lines from the four intersection points to any other point on first conic, and the second by finding the cross ratio of the lines from the same four intersection points to any other point on the second conic.
Therefore, let {ā k } 4 k=1 = A i ∩A j be the four intersection points of the conics in the image plane and letp be an arbitrary point on conic A i . Since the points are written in homogeneous coordinates, we may form the line joiningp andā k as w k =p ×ā k . Thus, the cross ratio Cr (w 1 , w 2 , w 3 , w 4 ) for ellipse A i is a projective invariant, and is the same for any choice of p on the ellipse. A second invariant may be computed Fig. 17 Chasles' theorem states that the cross-ratio formed by the family of red lines is the same as the cross-ratio of the family of blue lines in exactly the same way by picking an arbitrary point on ellipse A j . This cross-ratio may also be rewritten in terms of the eigenvalues of the matrix A −1 i A j [87]. Explicit computation of these cross ratios is not efficient, so an alternate technique is more appropriate-though the two alternate invariants are rational functions of the cross ratios.
A more conveniently computable form of the invariants for a pair of coplanar conics may be found by direct analysis of the matrices A i and A j describing the image conics. To see this, consider the pencil of conics λA i + μA j parameterized by the scalars λ and μ. Recall now that 3D crater C i projects to image ellipse A i according to the homography from Eq. 58. Since the craters are assumed coplanar their appearance in an image is related by a common homography, Taking the determinant yields which may be expanded as Thus we find that The coefficients Θ 1 , Θ 2 , Θ 3 , Θ 4 are only semi-invariants since their numerical value changes with |H C | and the arbitrary choices of s i and s j . These may easily be turned into rational invariants by considering appropriate ratios. To see this, first express the coefficients in terms of A i and A j , Likewise, we find that As shown in [128], the simplest pair of independent rational invariants are which we quickly verify produce a pair of non-trivial scalars that are not dependent on H C , s i , or s j (hence, they are the two fully generic invariants for a pair of conics): Since the left-hand side of these expressions are independent of H C , the same value is computed for any arbitrary projection of the conics C i and C j . These results are in agreement with [115], which also follows [128].
Following the convention of [41], we now observe that the equations simplify considerably by choosing the scale such that |A i | = |A j | = 1, making the coefficients Θ 1 and Θ 4 trivial. Using this simplification, substitute (108) and (111) into Eq. 106 to find which, in turn, leads to Therefore, when conics are pre-scaled to |A i | = |A j | = 1, the two unique invariants are simply We note that [41] follows a completely different approach to finding I ij and I ji , and even more geometrically-inspired derivations may be found in [96]. The specific approach shown here is chosen because of its extensibility to finding invariants of three or more coplanar conics. These two invariants for a pair of coplanar conics were used by Cheng, et al., for crater identification in [18,20]. We find, however, that better performance is often achieved by considering a triad of craters.

Triads of Coplanar Conics
It is established in Section 5.2 that there exist 5d − 8 projective invariants for a dtuple of coplanar conics (d ≥ 2). Thus, for a triad of coplanar conics, there exists 15 − 8 = 7 projective invariants. This fact has been known for some time [53,116].
The procedure for finding these seven invariants follows the same framework used for a pair of coplanar conics. Therefore, consider the determinant of a net of three conics [128]. As observed in [116], this evaluates to We immediately see that the coefficients of the first nine terms correspond to the pair-wise combinations of the three conics, with results that follow directly from the two-conic case and that Note that Θ 10 is the only coefficient that simultaneously depends on all three of the conics, thus making it the only term unique to a triad (and not to a pair). This may be computed as where the reader is briefly reminded that A * is the adjugate of A. Computation of this result (which, interestingly, does not appear in [116], [53], or any other reference discussing triads of coplanar conics) is tedious but straightforward, and is left as an exercise to the reader. The ten coefficients of Eq. 122 may be used to define the seven unique invariants. As with the pair of coplanar conics, the simplest approach is to choose the arbitrary scale of the matrices A i , A j , and A k such that |A i | = |A j | = |A k | = 1. Thus the coefficients Θ 1 , Θ 4 , and Θ 7 become trivial, and the remaining seven coefficients become the unique invariants. That is, with |A i | = |A j | = |A k | = 1, we find that Therefore, if three coplanar craters C i , C j , and C k project into three image conics A i , A j , and A k (described by the matrices A i , A j , and A k ), then the seven scalar values in Eqs. 128 to 131 will remain exactly the same under perspective projection regardless of the camera position and attitude.

Creating a Crater Pattern Descriptor from Invariants
Having identified projective invariants that may be computed from the observed crater rim geometry in an image (three for a triad of non-coplanar conics on a quadric surface and seven for a triad of coplanar conics), we wish to use this information to construct a pattern descriptor. This descriptor must have a structure that can be consistently reproduced to allow for matching against a static index, which requires some care since the invariants from Section 6 depend on the ordering of the observed craters. We note that the sensitivity of invariants to ordering is not inherently a bad thing (after all, we must specifically match each individual crater to the catalog)-though it does result in a few different methods for building descriptors. The index is built using descriptors for crater triads. We must choose, therefore, whether we want to impose ordering of the three craters within the triad before or after matching to the index. If we choose the former, then the invariants developed before may be used directly. If we choose the latter, we must either sort the projective invariants from Section 6 into a prescribed order or transform them into projec-tive and permutation (p 2 ) invariants. These conventions each have their advantages, which we now discuss.

Crater Pattern Descriptors with Projective Invariants
Suppose we observe a triad of craters that we want to match to a precomputed index in an order-dependent fashion. Since the crater pattern must be observed from above (i.e., the camera must be looking downwards from above the lunar surface because it cannot be inside the Moon), there are three possible orderings of the craters. If we choose to arrange these in a clockwise fashion within the image (see Fig. 18) then the three possibilities are (1, 2, 3), (3, 1, 2), and (2, 3, 1). For a match to be successful, the observed descriptor must match the index descriptor. Therefore, if the descriptor entries are order dependent, then a match only occurs when the observation ordering matches index ordering.
There are at least two options for matching to an index with order-dependent descriptors. First, if each crater triad has a single index entry, the index must be searched three times (once with each possible observation ordering). Alternatively, the index could include three entries for each crater triad (one for each ordering) and the larger index is searched only once. The first method (one index per triad) performs best in the presence of measurement noise.
The pattern descriptors in this scenario are simply a concatenation of the invariants from Section 6. Thus, for a triad of conics on a quadric surface, we build a threeelement descriptor: where the individual elements are from Eqs. 94 to 96. Likewise, for a triad of conics on a plane, we build a seven-element descriptor: I T ij k = I ij I jk I ki I ji I kj I ik I ij k (133) where the individual elements are from Eqs. 128 to 131.

Crater Pattern Descriptors with Sorted Projective Invariants
As before, suppose we observe a triad of craters. Now, instead of having an orderdependent descriptor, suppose that we want to match to a precomputed index in Fig. 18 Crater triad labels in an image are assigned in a clockwise fashion. There are three possible permutations of these labels, which may be interpreted as a cyclic permutation of the labels an order-independent fashion. There are a variety of approaches one might use to achieve this objective. We will see, however, that some care is required to arrive at a robust solution.
A very common technique for achieving permutation invariance is to always permute the order of the descriptor elements such the smallest (or largest) element appears first. For example, one possible descriptor is J T ij k = [J i , J j , J k ], where J i ≤ J j ≤ J k is a rearrangement of (J i , J j , J k ), i.e., we sort the numbers J i , J j , J k from small to large and let the result be the descriptor. This arrangement was employed by Hanak [49,50] to define the clockwise/counterclockwise sense of a crater pattern by ordering the legs of the crater triangle from shortest to longest. Similarly, for their lost-in-space algorithm, Maass, et al., [86] presort the entries in their index of crater triangles (which consists of two triangle interior angles) such that the smallest interior angle is first. Other examples of this idea abound within the space navigation community, ranging from crater pattern identification to star pattern identification.
Returning to the problem at hand, we know from our labeling convention that the crater pattern indices may only undergo a cyclic permutation. Thus, we may cyclically permute the entries until smallest valued invariant is in the first position-that is, where (J i , J j , J k ) is a cyclic permutation of (J i , J j , J k ) with J i = min(J i , J j , J k ).
Thus, for a triad of conics on a quadric surface, we build a three-element descriptor: where we label observed crater i such that J i = min(J i , J j , J k ) and the rest of the primed invariants follow by cyclic permutation. The individual elements are from Eqs. 94 to 96. Likewise, for a triad of conics on a plane, we build a seven-element descriptor: where we label observed crater i such that I ij = min(I ij , I jk , I ki ) and the rest of the primed invariants follow by cyclic permutation of the crater labels i, j, k (and, to be clear, not by cyclic permutation of the elements of the descriptor vector). The individual elements are from Eqs. 128 to 131. The usual argument for anchoring the descriptor element ordering with the smallest (or largest) element is that it is simple. We find in practice, however, that this ordering scheme is not optimally descriptive and is not very robust to noise. Indeed, it is not uncommon for noise to cause two similarly-valued invariants to switch position. If this happens to the value used to anchor the descriptor permutation then it renders the entire descriptor useless-even though the invariants themselves still describe the crater pattern well. Thus, we cycle through the sorted projective invariant descriptors exactly as we do for the unsorted case (Section 7.1). However, the presorting works most of the time, and greatly reduces the trials required (on average) to find a match.
The main problem with this approach is how sorting populates the k-d search space. By sorting the descriptor to have the smallest element first, we preferentially populate the search space for smaller values along this dimension. This doesn't matter so long as the measurement noise is very small-but, when supplied with noisy measurements, this reduces the likelihood that the correct descriptor will correspond with the first nearest neighbor. The end result is that an index of sorted descriptors will yield fewer matches (as compared with Section 7.1) when the data is noisy. Thus, while it may provide speed-up for low noise cases, both matching performance and speed quickly degrade as noise increases.

Crater Pattern Descriptors with Projective and Permutation (p 2 ) Invariants
Our objective is to construct a descriptor that is invariant to the order in which the craters are observed. It should be no surprise, therefore, that invariant theory is (once again) a useful mathematical framework for achieving this objective. There is some precedence for this approach using the so-called projective and permutation (p 2 ) invariants. The p 2 invariants are discussed for a set of five points in [76] and for a pair of conics in [87]. An attempt was made to apply this to the problem of crater identification in [110], though this approach suffered from a number of mistakes (see Section 3.3). Within the context of this past work, we now introduce the first complete set of p 2 invariants for a triad of conics on a quadric surface or on a plane (e.g., p 2 invariants for a triad crater rims on the Moon).
We briefly remind the reader that the p 2 invariants introduced here are algebraic functions of the projective invariants from Section 6. Thus, the p 2 invariants do not represent any new (independent) information. They simply express the information we already have in a different way.

Non-Coplanar Crater Patterns (Three-Element Descriptor)
We first develop the p 2 invariants for a triad of conics lying on a quadric surface. We know there exists exactly three algebraically independent invariants for this case, and we are able to compute three such invariants J i , J j , J k directly from an image using Eqs. 94 to 96. Our goal is to compute rational functions of J i , J j , J k that are invariant to a cyclic permutation of the crater labels.
Consider, therefore, the triple of numbers {x, y, z}. The standard generators of the ring of polynomial invariants for the symmetric group, S 3 , of all coordinate permutations are the elementary symmetric functions (see [103,Chapter 4]) which are clearly invariant to permutations of {x, y, z}. To obtain the generators of the invariant ring for the cyclic group Z/3Z, which we may also view as the alternating group A 3 , we require one more invariant: (see for example [103,Example 4.24]). We observe that Δ is invariant under cyclic permutations, but that odd permutations change its sign. While these polynomials could certainly be used to construct p 2 invariants and build a pattern descriptor, we note that the polynomials are of different degree which often leads to undesirable numerical properties for indexing. Thus, we instead look for three rational functions, where the difference in the degrees of the numerator and denominator are the same for all three functions.
To do this, we begin by defining a polynomial map F (x, y, z) = (F 1 (x, y, z), F 2 (x, y, z), F 3 (x, y, z)) (140) where F 1 , F 2 , F 3 are continuous rational functions with the property that F (x, y, z) = F (x , y , z ) if and only if (x , y , z ) is a cyclic permutation of (x, y, z). Let the cyclic group Z/3Z act on R 3 by permuting the coordinates. Then we desire F 1 , F 2 , F 3 to be invariant functions on R 3 . Whenever a finite group acts by a linear transformation, there exists a particular coordinate change where the action of the group becomes diagonal. That is, in some coordinate system, the matrix describing the action of each group element is a diagonal matrix. This follows from two facts from the representation theory of finite groups (see [42,Chapter 1]). First, every representation of a finite group is a direct sum of irreducible representations. Second, every irreducible representation of a finite abelian group is 1-dimensional. For the specific problem under consideration here (the action of Z/3Z on R 3 ), we can make this very explicit: namely, the change of coordinates making the action diagonal is a discrete Fourier transform (DFT) (which involves complex numbers). To make these ideas explicit, let a, b, c be where ζ = e 2πi/3 = − −1+ √ 3i 2 . Then the generator σ of Z/3Z acts by the diagonal matrix ⎡ We define F 1 , F 2 , and F 3 as rational functions of a, b, c, The denominator of F 2 and F 3 is equal to 1 2 and is equal to 0 if and only if x = y = z. We can also verify that This implies that if (x, y, z) is a sequence that converges to a point (t, t, t), then F 2 (x, y, z) and F 3 (x, y, z) converge to 0. Consequently, we can extend F 2 and F 3 to continuous functions on all of R 3 by defining them to be 0 whenever x = y = z. F (x , y , z ) if and only if (x , y , z ) is a cyclic permutation of (x, y, z) for all (x, y, z), (x , y , z ) ∈ R 3 .
Proof Suppose now that F (x, y, z) = F (x , y , z ). Let a, b, c be as given in Eqs. 144-146 and define similarly a = x + y + z , etc. To begin, observe that if F 2 (x, y, z) = F 2 (x , y , z ) = F 3 (x, y, z) = F 3 (x , y , z ) = 0 then we must have x = y = z and x = y = z and from F 1 (x, y, z) = 3x = 3x = F 1 (x , y , z ) it follows that (x, y, z) = (x , y , z ). Now suppose that F 2 (x, y, z) 2 + F 3 (x, y, z) 2 (and, therefore, This implies that Therefore b = ζ s b for some s ∈ {0, 1, 2}. Taking complex conjugation, we get that c = ζ −s c. We conclude that (a, b, c) and (a , b , c ) lie in the same orbit Z/3Z, which means that (x , y , z ) is a permutation of (x, y, z).
One can also verify that the field of rational invariants R(x, y, z) Z/3Z is generated by F 1 , F 2 , F 3 .

Coplanar Crater Patterns (Seven-Element Descriptor)
The development of p 2 invariants for a triad of conics lying on a plane follows similar arguments as for conics on a quadric, though it requires more care due to additional relationships between the seven elements in the descriptor. Therefore, recall from Eqs. 128-131 that the seven invariants consists of three invariant pairs, (I ij , I ji ), (I jk , I kj ), (I ki , I ik ), and a 7th invariant for the entire triad I ij k . Because they're constructed from a common set of crater observations, the two invariants within any given pair always undergo a common permutation. We also observe that I ij k is already symmetric and invariant to the ellipse label order.
Simply by applying the results of the previous section, we can find six Z/3Zinvariant functions from F (I ij , I jk , I ki ) and F (I ji , I kj , I ik ). The descriptiveness of these six invariants, however, is not optimal because members of a common pair are not forced to undergo the same permutation. Therefore, we introduce two new invariants.
Again making use of a DFT, define a , b , c as for = 1, 2. Note that c = b where x , y , z are real. A generator of Z/3Z acts on (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ) by the same cyclic permutation of coordinates, and it acts on b 1 , b 2 with a scalar ζ and on c 1 , c 2 with a scalar ζ 2 . We see, therefore, that b 1 c 2 is invariant. We define G 1 (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) and G 2 (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) be the real and imaginary part of b 1 c 2 , which we compute as We may isolate the real part to obtain G 1 and the imaginary part to obtain G 2
It follows, therefore, that there are many approaches to choose a descriptor for a set of three coplanar craters. We summarize three obvious choices, with the third being the one we usually recommend.
The most obvious approach is to form a seven-element descriptor as

[F (I ij , I jk , I ki ), F (I ji , I kj , I ik ), I ij k ]
As was noted above, the descriptiveness of such a scheme is not optimal, because this descriptor remains the same when the values of I ij , I jk , I ki cyclically rotate while the values of I ji , I kj , I ik stay the same. This means that 3 different configurations might have the same descriptor. A second option is to form a nine-element descriptor

[F (I ij , I jk , I ki ), F (I ji , I kj , I ik ), G(I ij , I jk , I ki , I ji , I kj , I ik ), I ij k ]
This descriptor is more descriptive but uses more space. Here, we have replaced the invariant G with G, which we define as This is an essential step, since G scales quadratically with (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) while F scales linearly. Using G instead of G results in a poorly scaled descriptor that complicates (and sometimes precludes) nearest neighbor searches with efficient data structures. Conversely, the function G scales linearly. We note that G is not a rational function because of the 4-th root, but it does extend to a continuous function on R 6 by defining the function to be 0 whenever x 1 = y 1 = z 1 or x 2 = y 2 = z 2 . A third option (and the one we suggest) is to use the seven-element descriptor For generic configurations, this descriptor is as good as the nine-element descriptor. The subfield of L = C(x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) generated by F (x 1 , y 1 , z 1 ), = 1, 2, 3, F 1 (x 2 , y 2 , z 2 ) and G (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ), for = 1, 2 is equal to K = C(a 1 , a 2 , b 3 1 , b 1 c 1 , b 1 c 2 , c 1 b 2 ). We observe that L = K(b 1 ) and b 1 has degree three over K. So the degree of the field extension L/K is at most three, K is contained in the fixed field L Z/3Z and the degree of the extension L/L Z/3Z is three. This implies that K = L Z/3Z . In other words, every Z/3Z-invariant rational function in I ij , I jk , I ki , I ji , I kj , I ik is a rational function in the first 6 elements of the descriptor.
However, in the degenerate case I ij = I jk = I ki = t for some fixed t we have G(I ij , I jk , I ki , I ji , I kj , I ik ) = 0. In that case, the descriptor cannot distinguish between different configurations of the same pattern (e.g., a pattern and it's mirror). This is to be expected, since the degenerate case corresponds to equally sized craters around an equilateral triangle. Here, while there may be enough information to uniquely identify the triad (i.e., match the triad descriptor to the database), there is not sufficient information to disambiguate the specific crater labels with the projected crater rims alone. No special action is required, since this ambiguity is naturally handled in the pattern verification process (Section 9).

Remarks on Choosing a Descriptor Convention
In the sections above we introduce descriptors based on the projective invariants, sorted projective invariants, and on the p 2 invariants. Which one is best is often application dependent.
Descriptors built directly on the projective invariants (Section 7.1) have better matching performance when presented with noisy data (see left-hand plot of Fig. 19), though they require the index be searched three times per observed pattern (one for each possible permutation).
Conversely, the descriptors built with the sorted invariants (Section 7.2) or p 2 invariants (Section 7.3) only require one index search per observed pattern, but are more sensitive to measurement noise. As a consequence the sorted invariant and p 2 invariant descriptors are faster for low noise situations, but slower for high noise situations (see right-hand plot of Fig. 19). This may seem counter-intuitive. The explanation, however, is straightforward. As measurement noise increases, the sorted invariant and p 2 invariant descriptors must attempt more triads (on average, as compared to the projective invariant descriptor from Section 7.1) before finding a correct nearest neighbor match. Thus, with large amounts of measurement noise, the projective invariant descriptor from Section 7.1 tends to find a match sooner despite needing three index searches per triad-and finding a match sooner results in a lower (faster) run-time.

Building a Global Crater Index
Global lunar crater identification requires multi-scale indexing and careful catalog curation. When close to the Moon, a spacecraft sees small craters that are nearly coplanar. Conversely, when far from the Moon, a spacecraft sees larger craters distributed over the lunar sphere (not coplanar). This immediately suggests at least two crater indexes be built-one for small patterns of nearly coplanar craters and another for larger patterns of non-coplanar craters. In practice, we find that the most efficient real-time performance is achieved by creating a single index for the planar case and Fig. 19 A comparison of crater matching performance and execution speed for 1,000 random images of the Moon. These images correspond to a spacecraft at 600 km altitude placed randomly around the Moon (uniform distribution over the lunar sphere). Matching is performed using the three-element (non-coplanar) crater triad descriptor from Section 6.1.3. The left-hand plot shows the fraction of cases with a successful match. The balance of cases returned no match, and there were no cases where an incorrect match was returned. The right-hand plot shows execution time for the entire matching pipeline (conics as input, verified match as output) using the three different descriptors. These times correspond to non-optimized prototype code and should be taken to represent relative (rather than absolute) execution time a hierarchy of indexes for the non-coplanar case. These indexes may be combined into a single large index or kept as separate indexes. We find the latter to often be the better choice.
The lunar crater database from Ref. [119] contains 1.3M craters over about 1-2 km in diameter. It is immediately obvious that consideration of all the 1.3M 3 = 3.66 × 10 17 combinations of possible crater triads is not reasonable. This further motivates construction of a hierarchy of indexes with different scales to better manage the combinations of craters that could plausibly be observed at the same time. A visualization of the local and global crater indexes used in this work is shown in Fig. 20, where the index spaces are superimposed on the lunar crater density data from Ref. [119].
We propose the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) 4 framework [45] be used to subdivide the lunar surface into equal area regions, which we refer to as surface pixels. HEALPix was developed for subdividing the celestial sphere in support of science objectives for the Wilkinson Microwave Anisotropy Probe (WMAP) mission [10], and has since found widespread use for analysis of data from other NASA and ESA missions [e.g., Cosmic Background Explorer (COBE), Planck]. It is also used for managing the creation of star quadrilaterals in the present state-of-the-art for star identification, calibration, and alignment of astrometric images [71]. We believe our work to be the first application of HEALPix to the management of lunar surface features.
All indexes are built using the same fundamental approach, but at different scales. The procedure is as follows: First, the lunar surface is tiled into N pix = 12(2 2k ) surface pixels of equal area. Next, a list is constructed for each pixel containing the catalog entries for craters within a specified size range (e.g., minimum/maximum diameter) and whose catalog fit was constructed using at least 90% of the rim's circumference.
Second, given a list of usable craters in each HEALPix surface pixel, we loop through all the pixels to create crater triads. At each pixel, we consider craters from the 3 × 3 HEALPix grid centered about the reference pixel. From these 9 surface pixels, all possible triads are formed where (1) craters do not intersect one another and (2) the triad center lies within the reference (center) surface pixel. Valid crater triads are arranged in a clockwise order and the scale-appropriate projective invariant descriptors are computed. To better illustrate this, Fig. 21 shows an example 3 × 3 HEALPix grid (region of support for surface pixel 6318) overlayed on an example Metric Camera image from Apollo 17. For further intuition, a global visualization of the HEALPix grid and the crater density (i.e., number of craters per HEALPix surface pixel) is shown in Fig. 22.
Third, once the triads belonging to all surface pixels have been computed, the results are stored in an efficiently searchable data structure. The authors have found excellent performance using either a k-d tree [11] or n-d k-vector [6], though other reasonable choices exist. We briefly remark that the 1-d k-vector [93,94] has seen  [73,121] extensive use for star pattern matching in space applications [95,134], and the n-d kvector has been proposed for space applications as well [74]. The specific choice of data structure is not of primary concern in the present analysis so long as a reasonable selection is made.

Indexing Local Crater Patterns with Coplanar Invariants
To construct a global index of local-scale (small) crater patterns, we choose a HEALPix resolution of k = 5 for craters of diameter 4-30 km. This yields N pix = 12(2 10 ) = 12, 288 surface pixels, with each surface pixel having a surface area of approximately 3, 086 km 2 . There are 20,737 craters with catalog parameters supported by > 90% of the crater rim circumference and having a diameter between 4 km and 30 km. Rather than producing 20,737 3 = 1.49 × 10 12 combinations, the HEALPix grouping strategy only produces the 4.8M crater triads that are nearly coplanar.
The local crater patterns are assumed coplanar and craters are allowed to follow an arbitrary elliptical shape. Thus, we construct a seven-element descriptor using an appropriate method from Section 7.

Indexing Global Crater Patterns with Non-Coplanar Invariants
For crater patterns that are regional or global in extent, we find it helpful to construct at least two indexes-though specific mission needs may require more or less. The extension to a single global index or many global indexes is trivial. The Moon is very nearly a sphere at the global level. Therefore, because the noncoplanar invariants require each pair of ellipses to lie on a common quadric surface, the craters must be nearly circular. The global lunar crater databases constructed here are built only from craters having an ellipticity of a/b ≤ 1.1. Since the global patterns are assumed to lie on a nondegenerate quadric surface, we construct a three-element descriptor using an appropriate method from Section 7.
The first global crater index has a HEALPix resolution of k = 3 for craters of diameter 25-125 km. This yields N pix = 12(2 6 ) = 768 surface pixels, with each surface pixel having a surface area of approximately 49, 390 km 2 . The index is good for matching crater patterns at the regional level. There are 904 nearly circular craters (a/b ≤ 1.1) with catalog parameters supported by > 90% of the crater rim circumference and having a diameter between 25 km and 125 km. Rather than producing 904 3 = 1.23 × 10 8 combinations, the HEALPix grouping strategy only produces the 140,929 crater triads. The example crater triad in Fig. 16 (see Section 6.1) was recognized using this index.
The second global crater index has a HEALPix resolution of k = 1 for craters having a diameter over 100 km. This yields N pix = 12(2 2 ) = 48 surface pixels, with each surface pixel having a surface area of approximately 790, 300 km 2 . The index is good for matching crater patterns at the global level, when nearly an entire hemisphere is visible in an image. There are 31 nearly circular craters (a/b ≤ 1.1) with catalog parameters supported by > 90% of the crater rim circumference and having a diameter over 100 km. Rather than producing 31 3 = 4, 495 combinations, the HEALPix grouping strategy only produces the 707 crater triads.

Remarks on the Utility of Crater Index Hierarchies
The primary purpose of having more than one index is to apply coplanar invariants at the local scale (allowing for arbitrary elliptical crater shape) and non-coplanar invariants at the regional/global scale (where curvature of the Moon makes craters lie in substantially different planes).
In practice, the orbital regime is often known ahead of time (e.g., LLO, cislunar), allowing us to only search crater patterns in the index of appropriate scale. If, however, we are truly "lost-in-space" and do not know the orbital regime a priori, we have found that all indexes may be queried and only the correct one will produce a match. Thus, the hierarchy of indexes allows us to exploit mission-specific information when it is available, but does not necessarily require this information.
While this work builds a truly global index for both local crater patterns and global crater patterns, it may sometimes be desirable to only include craters visible along a reference trajectory.

Matching Observed Crater Patterns to a Pre-Built Index
When supplied an image of the lunar surface, the objective is to recognize a pattern of observed craters using the image conics corresponding to the projection of the crater rims.
This is accomplished using a straightforward matching and and verification procedure. First, a crater detection algorithm (CDA) produces an ellipse fit to a crater rims in an image. Second, we consider triads of these image conics, cycling through combinations using the Enhanced Pattern Shifting (EPS) method [5]. For each crater triad, we compute the appropriate descriptor (see Section 7) and query the corresponding pre-built index to find potential matches (e.g., best N ≥ 1 matches). For each crater match hypothesis, we compute the unknown camera position (see Section 9.1). After first checking to ensure the position estimate is not inside the Moon (as would happen if we accidentally match to a pattern on the wrong side of the Moon), we then reproject the expected crater rim locations in the image using Eq. 47. The match is verified by comparing the observed and expected crater rims using a new distance metric (see Section 9.2). Previous methods only compared crater centers and thus required 4-5 correct correspondences to verify a match. Comparing the craters rims is more informative, allowing for pattern match verification with only three craters (thus allowing crater identification above more sparsely cratered terrain).
If a pattern match hypothesis is verified by this procedure, we declare this the solution and terminate the search. If a pattern match hypothesis is not verified we consider another hypothesis for that triad or move to the next triad in the EPS sequence. If we reach the end of the EPS sequence with no hypothesis being verified, we declare that no match is possible and terminate the search.
The subsections that follow describe the mathematics for pose estimation from corresponding non-coplanar conics and a crater rim distance metric. Both of these algorithms are novel.

Pose from Non-Coplanar Conics in Correspondence
The usual approach for computing pose from matched craters only uses the coordinates of the crater center. Here, it is essential to remember that the center of an image ellipse does not generally produce a line-of-sight vector to the center of a 3D ellipse (or circle). Some past crater identification pipelines consider this effect (e.g., [110]), but many do not. Regardless, there are a variety of algorithms one may use to compute pose from corresponding 2D image and 3D model points [4,77]. Unfortunately, however, using only the crater center points neglects the substantial amount of valuable navigation information contained within the shape of the projected crater rims. We aim to address this problem here.
Were the 3D conics strictly coplanar, there are a number of algorithms one could use to estimate pose from the projected image conics [67,135,141]. However, while we assume local crater patterns are coplanar for index building/matching, there is no reason to assume coplanarity for pose estimation or match verification. Furthermore, we require the ability to solve this problem for both local and global patternsthus requiring a method for pose estimation from the projection of non-coplanar conics.
The literature discussing pose estimation from non-coplanar conics is limited. The most common approach is to look at a pair of ellipses in an image and construct a single-parameter family of possible poses from one of these two ellipses. This one parameter family is then numerically searched to find the best possible agreement with the second ellipse. This procedure was suggested in [85], which is generally cited as the solution to this problem.
The approach of [85], however, is not appropriate for spacecraft navigation as it (1) does not fully use the information content of every conic and (2) it only deals with a conic pair. Instead, we search for a solution that uses all the information from a d-tuple (usually a triad) of image conics to solve for camera pose in the least squares sense.
Since the application at hand is lunar crater identification, we may simplify the problem by assuming the relative attitude is known. We consider this to be a reasonable assumption in practice since we have excellent knowledge of both the Moon's attitude (from ephemeris data; e.g., SPICE kernels [1,2]) and the spacecraft attitude (from star trackers [81,82]). There are two observations that can now be made. First, we can think of no plausible failure mode (where recovery is still possible) in which the spacecraft has no knowledge of time (necessary for finding lunar attitude from SPICE kernels) or of inertial attitude. Second, we observe that this still formally qualifies as a lost-in-space problem since it assumes no knowledge of the spacecraft translational states (position or velocity). Moreover, the "known" attitude comes from the star tracker and there are many well-established lost-in-space star identification algorithms [118,134]. Regardless, we presume the attitude transformation matrix T M C is known, such that the pose estimation problem only requires a solution for the selenographic camera location r M .
Therefore, let us begin with consideration of the projection of a single conic as described by Eq. 58. Now, substitute for H C i from Eq. 50 and expand to find where Recall here that A i describes the measured image conic, C i describes the 3D crater conic in the ith crater's ENU frame, T M C is presumed known, and k T = [0 0 1]. Thus, the only unknown in Eq. 150 is r M . Now, recalling (34) for H M i , rearrange the left-hand side into matrix form where S is from Eq. 35. It is immediately evident that the upper-left 2 × 2 submatrix is independent of r M , the lower-right element is quadratic in r M , and the remaining terms are linear in r M . Therefore, we first find the unknown scalar s i using only the upper-left 2 × 2 submatrix. Using the Frobenius norm, computeŝ i aŝ which has the least squares solution With the scale s i known, we can now take the upper-right 2 × 1 submatrix to form the linear system for all of the observed conics This may be solved for the selenographic camera location r M in the least squares sense.

Comparing Two Image Conics
In order to verify a crater match hypothesis, we require a measure of the distance between two conics. For crater pattern verification, this distance of interest is usually between the image ellipse we expect from the projection of a crater's rim (A i ) and the image ellipse we measure of the same crater's rim (Ã i ). For the moment, however, we will briefly discuss how to compute the distance between two arbitrary ellipses: Given two ellipses in an image, A i and A j , we seek a scalar distance metric d(A i , A j ) that satisfies the three usual axioms for a distance metric [32] 1. Minimality: d(A i , A j ) = 0 iff A i = A j . That is, the distance between an ellipse and itself is zero.
as well as a fourth axiom unique to this application 4. Similarity Invariance: is a similarity transformation. That is, the distance between A i to A j should not change if the two ellipses undergo a common translation, rotation, or scaling in the image (i.e., undergo a common similarity transformation). See Fig. 23.
After reviewing the literature, we found numerous approaches for comparing ellipses, but all fail to meet one (or more) of the above four axioms. Thus, after a brief review of existing techniques, a novel method is proposed. There exist various ad hoc methods based on the explicit ellipse parameters (a, b, x c , y c , and ψ) [30,113] or by measuring the distance only at one specific point (e.g., via the Hausdorff distance [132,136]), though none of these are similarity invariants. Moreover, these methods have changing geometric meaning as the two ellipses change shape and relative orientation. While not a violation of one of the axioms, it is not a desirable attribute.
Another approach is to apply a common normalization to the implicit representation of both conics (e.g., det(A i ) = 1) and then compute the Frobenius norm of their Fig. 23 Examples of an ellipse pair undergoing a simiarity transformation. The left-most ellipse pair may be translated (S 1 ); translated and rotated (S 2 ); or translated, rotated, and scaled (S 3 ). Under the similarity invariance axiom for an ellipse pair distance metric, all four of these ellipse pairs should produce the same value for distance since they all have the same relative geometry. Observe that all four of these ellipse pairs have the same Jaccard distance and Gaussian angle. The area of ellipse intersection (i.e., ellipse overlap) is shaded in light blue difference [59]. There are various normalizations one could choose for the matrix A, as summarized in [66]. Nevertheless, such a comparison metric is given by though the specific geometric meaning of this is not clear. This metric also fails to meet the similarity invariance axiom.
We also considered computing the distance as the line integral about one ellipse of the Euclidean distance to the other ellipse. While this admits an elegant and efficient solution with many interesting properties, it fails to meet both the symmetry and similarity invariance axioms. To meet the symmetry requirement, it would be possible to compute this metric in both directions and sum the result, though this is expensive. Regardless, a distance computed in this fashion will still fail to meet the similarity invariance axiom. The method was thus abandoned for the present application.
The deficiencies of the above methods motivate the need for another metric that satisfies all four distance axioms, has a geometrically consistent meaning, and possesses well-understood statistics. There are at least two such distance metrics: one using the Jaccard distance and another by interpreting the ellipse as a multivariate Gaussian.

Ellipse Pair Distance Metric with the Jaccard Index
The Jaccard index is the ratio of an intersection to a union. In this case, we take it to be the ratio of the areas of the ellipse intersection to the ellipse union, The Jaccard distance, which is well-known to satisfy all three of the classical distance metric axioms [79], is computed as It is straightforward to see that this metric is also satisfies our similarity invariance axiom.
The only difficulty with the Jaccard distance is the quick computation of the intersection area. Although it is possible to analytically compute the area of ellipse overlap [60], doing so requires the separate consideration of nine different relative ellipse configurations. Such algorithm branching is generally undesirable. Fortunately, we recall that the ellipses in question represent regions of substantial extent in a digital image. Thus, the Jaccard distance may be approximated as the ratio of the number of pixels inside both ellipses to the number of pixels inside either ellipse. We may quickly count the pixels inside A i by finding those pixels that satisfy the inequalitȳ u T A iū < 0, which is a simple and parallelizable computation. Once we count the pixels in this way, computing the Jaccard distance is trivial.

Ellipse Pair Distance Metric with the Gaussian Angle
It is possible to interpret the 2D ellipse A i as a bivariate Gaussian probability density function. If that interpretation is used for both A i and A j , the distance between two ellipses may be computed using any number of measures for comparing Gaussian PDFs. The authors of [144] follow this interpretation and consider the use of the Kullback-Leibler Divergence (KLD) and variants of the Wasserstein distance. The KLD, however, does not satisfy either the symmetry or triangle inequality axioms. The Wasserstein distance meets the first three axioms, but is not invariant under scaling (hence, fails the similarity invariance axiom-though modifications to the classical metrics may alleviate this shortcoming).
We introduce a different distance metric also based on the interpretation of A i as a bivariate Gaussian. Therefore, begin by subdividing the conic locus matrix A i describing the image ellipse A i as and it follows that thatū T i A iūi = 0 and that With this particular scaling of A i , we immediately recognize the 2 × 2 submatrix Y i as the shape of the ellipse and y T i = [u c i , v c i ] as the image coordinates of the ellipse center. If we interpret the ellipse as the 1σ isofrequency contour of a bivariate Gaussian distribution (i.e., the so-called 1σ error ellipse), the ellipse A i corresponds to the distribution N(y i , Y −1 i ) with probability density Therefore, given two ellipses A i and A j , we may compute the distance between these two ellipses as the "angle" between p A i (u) and p A j (u). 5 Such an angle is computed as the inner product To compact integral notation in the equations that follow, we will no longer explicitly write the argument of integration. To proceed, we compute a distance metric as From inspection, this clearly satisfies the minimality and symmetry conditions. Satisfaction of the triangle inequality follows from the Cauchy-Schwarz inequality. We compute To find the numerator, we first write Now, considering the term inside the exponential, Thus, noting that c is a constant, we rewrite the integral as Consequently, We may therefore compute the distance metric d GA by substitution of Eqs. 165 and 174 into Eq. 164, In addition to satisfying the four required axioms, one of the great advantages of this particular distance metric is that it may be analytically computed from the parameters of the two image ellipses A i and A j .

A Fast Method for Assessing Ellipse Correspondence
We propose to use the Gaussian angle distance metric as the criteria for evaluating a potential crater correspondence. In the case where A j is a perturbed version of A i , we find that, to an excellent approximation, where The validity of this approximation is illustrated in Fig. 24 for both a nearly circular crater and a very elliptical crater. Since the distance metric is observed to follow a χ 2 4 distribution, we may construct a statistically-informed criterion for accepting a crater match hypothesis. This is accomplished by evaluating the quantile function for the χ 2 4 distribution. For example, the 99th percentile for the χ 2 4 distribution occurs at 13.277, and would lead to an acceptance criteria of The appropriate percentile for accepting or rejecting craters is application dependent and left as a design variable left to the analyst. Finally, when selecting the threshold for accepting crater matches, the user may wish to consider that very elliptical craters have heavier tails (as compared to a true χ 2 4 distribution) -thus increasing the likelihood of rejecting a valid crater assignment. This may be seen in the right-hand example of Fig. 24. Except for very elliptical craters, the deviation from the χ 2 4 distribution is modest and numerical experiments can help the analyst decide if this is important for a specific application.

Sensitivity to Crater Rim Fit and Viewing Geometry
We performed a Monte Carlo analysis to evaluate the sensitivity of the proposed crater matching methodology to differing levels of measurement noise and under different viewing geometries. Therefore, as an example, consider a spacecraft placed randomly (with a uniform distribution) around the lunar sphere, thus creating a situation where we randomly image with equal probability any part of the lunar globe. We assume a camera with a full field-of-view (FOV) of about 73.7 × 73.7 deg (same as the Apollo metric camera [38,143]) and with a 2, 200 × 2, 200 pixel focal plane array (creating a 4.8 MPixel image).
We now consider two parametric scenarios: (1) a nadir pointing camera with varying errors in crater rim localization and (2) a fixed level of error in crater rim localization and increasing off-nadir pointing angle. For the first scenario (nadir pointing images), let the error in the ellipse parameters a, b, u c , v c all vary by σ img and range from 0 to 3 pixel. Subpixel estimation of these parameters is not uncommon since we generally fit an ellipse to hundreds of pixel measurements around the crater circumference. Thus, for the second scenario (off-nadir pointing images), we assume an error of 0.5 pixel in the ellipse parameters.
Performance for each case is summarized with the root-sum-square (RSS) error in the estimated position of the camera. In all cases the estimated position errors appeared to be unbiased and are not explicitly recorded here for brevity.

Local Crater Matching Results
Suppose the randomly placed spacecraft has an altitude of 150 km above the lunar surface. At this altitude, the camera views primarily local crater patterns. In this case, results for the two Monte Carlo scenarios are shown in Tables 2 and 3.
From Table 2 we observe that matching performance is not appreciably affected by ellipse localization error until these errors reach about 2-3 pixels. Even then, the degradation in performance is gradual. The dropping percentage of correct matches occurs because the closest (nearest neighbor) index match to the computed invariants is no longer correct. We could delay the onset of this effect by considering additional neighbors beyond just the single best match, but this would increase run-time. The reader should also be warned that these results do not suggest that crater localization error always needs to be better than 2-3 pixels. What matters is not the absolute rim localization error, but rather the rim localization error as compared to the size of the A nadir-pointing camera is placed randomly around the lunar globe at an altitude of 150 km, with ellipse fit error varying from 0 to 3 pixel image craters (i.e., a 2 pixel error on a crater with a 10 pixel diameter is significant, while a 2 pixel error on a crater with a 500 pixel diameter is not significant). Test cases are counted in the "No Match Found" column when we see at least three craters but the algorithm produces no matches. This can happen for a number of reasons. First, the acceptance threshold for a crater match was set at the 99th percentile (as in Eq. 178), which would suggest that 1 − 0.99 3 ≈ 3% of the triads should fail to match due to measurement noise (this acceptance threshold can be adjusted based on the analyst's preferences). When there are more than three craters in an image, the likelihood of this happening becomes rather small (but is never zero). The second reason for no match is that no combination of observed craters corresponds to an entry in the index. Given the FOV of this particular camera, it is possible to see more than nine HEALPix surface pixels at one time-thus creating the possibility of observing three craters that do not correspond to a triad within the index.
From Table 3 we observe that off-nadir pointing has no meaningful effect (at least up to 30 deg) on our ability to recognize a crater pattern.   A nadir-pointing camera is placed randomly around the lunar globe at an altitude of 600 km, with ellipse fit error varying from 0 to 3 pixel

Global Crater Matching Results
Suppose the randomly placed spacecraft has an altitude of 600 km above the lunar surface. At this altitude, the camera views primarily regional/global crater patterns. In this case, results for the two Monte Carlo scenarios are shown in Tables 4 and 5. We observe the same basic trend for the global pattern as for the local patterns. Matching performance degrades as we approach localization error of around 2-3 pixels, which occurs at the roughly same level of noise only because the 600 km altitude choice. We once again observe very little affect on matching performance due to off-nadir pointing.
Note that the camera location is computed in exactly the same way for the local and global matching examples (both use the algorithm from Section 9.1). The global camera location errors are larger because the images are taken from a higher altitude. Regardless, errors are still less than 1 km for the local crater patterns and less than 5 km for the global crater patterns-even with rim localization errors much larger than should be expected in practice. A camera is placed randomly around the lunar globe at an altitude of 600 km, with an off-nadir pointing varying from 0 to 30 deg. Crater fit error of σ = 0.5 pixel on ellipse parameters

Local Matching on Clementine Images
We demonstrate the seven-element local matching descriptor on real images collected by the Clementine spacecraft's UV/Visible (UV/Vis) camera. During its 71 days in lunar orbit in 1994 [104], the Clementine mission collected millions of images of the lunar surface with a variety of instruments. One of these instruments was the UV/Visible camera, a 4.2 × 5.6 deg FOV camera with a 288 × 384 pixel CCD focal plane array [58,70]. Images of the lunar surface are available at a variety of ranges and viewpoints and serve as an excellent source of data for testing our proposed crater identification algorithms. Some example crater identification results are shown in Fig. 25. These results show good matching performance, both for nearly nadir pointing images (top two rows) and for images from oblique viewpoints (bottom row). For each example image, the best fit ellipse parameters (as manually obtained from the image) were perturbed by a Gaussian distribution with a standard deviation of σ img = 0.5 pixel. The noise causes a different crater pattern to be matched first for each run, resulting in different crater triads shown in Fig. 25. In each case, the craters with the white outline (and with centers connected by the thin yellow line) were correctly matched to the Robbins catalog [119] using the local crater index (HEALPix with k = 5) described in Section 8.

Global Matching on Synthetic Images on Reference Lunar Trajectory
We demonstrate the three-element global matching descriptor on synthetic images along a reference lunar trajectory. Fig. 26 Example synthetic images along a reference trajectory as the spacecraft moves away from the Moon (top to bottom). The craters used to compute invariants are outlined in white and a sampling of other craters (which are not used) are outlined in red. Monte Carlo results (1,000 trials) assume a 1-σ ellipse localization error of 1 pixel, leading to widening of invariant histograms as craters become smaller relative to the fixed error. Synthetic images of the Moon were produced using PANGU Planet Surface Simulation Software developed by the Space Technology Centre at the University of Dundee, Scotland [111] The synthetic images were generated using PANGU, and have a resolution of 2, 048×2, 592 pixels with a 30 deg horizontal FOV. Three example images are shown in Fig. 26. For each image, the parameters of reference crater ellipse (a, b, u c , and v c ) were then perturbed with uncorrelated Gaussian noise of σ img = 1.0 pixel. This was repeated 1,000 times for each image (with each trial having different noise applied to the crater rim fits) to understand the sensitivity of the invariants to measurement noise. Histograms of the three invariants are also shown in Fig. 26. Since the measurement noise is fixed (σ img = 1.0 pixel on a, b, u c , and v c ), we see the histograms widen as the spacecraft gets farther away from the Moon (top to bottom in figure). This occurs because the size of the pattern shrinks relative to the fixed magnitude of ellipse localization error-hence the fixed error causes a greater perturbation of the crater rim geometry.

Conclusion
The identification of craters in a digital image of the lunar surface is a critical capability, both for terrain relative navigation (TRN) and for the registration of scientific images. In this work, we provide the first comprehensive and mathematically rigorous approach for image-based crater identification using concepts from invariant theory. We first show that crater rims tend to be elliptical in shape, suggesting that they are well modeled as a conic (e.g., circle, ellipse) lying on the lunar surface. We then show that there are no projective invariants for conics lying on the surface of an arbitrarily shaped celestial body, but there are invariants for regularly shaped bodies like the Moon. Specifically, for a d-tuple of conics, we show that there are 3d − 6 algebraically independent projective invariants for conics lying on a nondegenerate quadric surface and there are 5d − 8 algebraically independent projective invariants for conics lying on a plane.
With the number of algebraically independent projective invariants known, we then develop a practical means of computing these invariants from just the apparent crater rims in a digital image (i.e., from the perspective projection of the 3D conics). These invariants may be used to form a pattern descriptor that is always the same, regardless of camera viewpoint. Thus, these descriptors may be precomputed for known crater patterns and stored in a searchable index. We find that descriptors for a triad of craters provide good matching results.
When given an image with observed crater rims, a match hypothesis for the crater rim pattern may be found by a simple nearest neighbor search of the index. With the resulting image-to-map crater correspondence hypothesis, we show how to estimate the location of the camera in the least-squares sense using the three crater rims directly (instead of just the crater center coordinates). Finally, to verify the crater match and pose hypothesis, we reproject the expected crater rim contours into the images and compare these with the measured crater rims. This is accomplished with a new distance metric that compares the entire ellipse fit (instead of just the crater center coordinates). By verifying with a distance metric considering the full crater rim contours, we are able to verify a pattern using fewer craters than methods using only center coordinates.
These techniques are demonstrated in a variety of numerical experiments-both on synthetic and real images. For a camera placed randomly around the Moon with no a priori position knowledge, successful matches generally exceed 90% of the test cases. The remaining cases return "no match," either because there were fewer than three craters in a particular image or because measurement noise corrupted the crater rim fits too much for unambiguous recognition. In no case was an incorrect match returned, a sign that the new distance metric can indeed verify match hypotheses with only three craters (instead of the usual ≥ 5 craters).
While this work represents a significant advance in methods for crater-based TRN, a great deal of work remains to be done. What follows are a few thoughts regarding the next steps: 1. While the repurposing of scientific crater catalogs for TRN is convenient, these catalogs were built with different objectives and are generally suboptimal for TRN applications. Thus, we suggest that the lunar exploration community would benefit greatly from a purpose-built TRN crater catalog. Such a catalog would be populated with only well-defined craters, with each having a shape and size consistent with what an automated CDA would produce. 2. We provide a complete set of independent projective invariants for crater patterns on the surface of the Moon in Section 6 (three invariants for a triad of conics on a quadric surface, seven invariants for a triad of conics on a plane). While there are no additional algebraically independent invariants, it is possible to construct other algebraically dependent invariants as functions of the ones provided in this work (e.g., the p 2 invariants in Section 7 are algebraic functions of the projective invariants from Section 6). There may exist other forms of these invariants that make lunar crater patterns more distinct or that have superior stability in the presence of measurement noise. Searching for such alternative formulations of the invariants from Section 6 is an interesting problem. 3. This work develops a global index of crater descriptors at three different scales using HEALPix. Most real missions likely do not need a global index-but, instead, just an index of craters along a reference trajectory. Although HEALPix may still be used to manage varying scales, there is a great deal of forward work in customizing the proof-of-concept outlined here to a mission-specific application. The best design for an actual in-flight index is not self evident and is the topic of ongoing work. 4. Not addressed in this work are crater detection algorithms (CDAs). Though many different CDAs exist, it is not clear which specific algorithms are best paired with the crater identification framework outlined here. Our crater identification algorithm is critically dependent on the localization of the best-fit ellipse to the crater rim, which is often not the metric used to design or evaluate CDAs. Thus, finding (or developing) an appropriate CDA is an obvious topic of follow-on work.
In conclusion, we find that crater-based TRN holds great promise for lost-in-space navigation near the Moon. We introduce a robust method for crater identification in this manuscript, though we note that additional work is required to achieve a flightcapable system. Crater-based TRN is a rich field for future study.