Authors: John A. Christian, Harm Derksen, Ryan Watkins
Categories: Original Article, Conics, Craters, HEALPix, Invariant theory, Lunar exploration, Pattern recognition, Projective geometry, Spacecraft navigation, Terrain relative navigation
Source: The Journal of the Astronautical Sciences
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.
Keywords: Spacecraft navigation, Invariant theory, Lunar exploration, Terrain relative navigation, Craters, Conics, Projective geometry, Pattern recognition, HEALPix
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-in-space 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 (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
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.
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 invariants1 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 pre-flight 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(logn) 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 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
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.
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 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 D90 is the crater diameter for a vertical impact and L is the impactor diameter (i.e., D90/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, D90/L∼60) and 30 deg (for brittle rock, D90/L∼4) [29, 39, 91].
Fig. 2 Geometry (side view) of impact crater creation, with the impact angle being defined as the angle above the local horizontal
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]
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 terrace-like structures. As impact features increase in size, they also increase in complexity. Craters larger than 300 km are usually classified as basins.
Fig. 4 Craters of different sizes and ages exhibit a variety of morphologies. (A) Lichtenburg B (33.3 N, 61.5 W), a 5 km simple bowl-shaped crater with a raised rim. (B) Aristarchus (23.7 N, 47.4 W), a 40 km complex crater with terraced walls and a single central peak. (C) Copernicus (9.62 N, 20.08 W), a 93 km complex crater with double central peaks, terraced walls, and a flat floor. (D) Riccius W (38.9 S, 25. 2 E), a 19 km degraded crater with an infilled cavity and much of the rim no longer clearly visible
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 ak (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.
Fig. 5 The Neukum Production Function (NPF) describes the cumulative number of new lunar craters per km per Gyr below a specified diameter. The red line shows the NPF evaluated at a diameter of 4 km
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.
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.
Of special note is the tremendous progress made in CDAs with machine learning over the last five years. Some notable algorithms DeepMoon2 [131], Python crater detection algorithm (PyCDA)3 [69], LunaNet [37], CraterIDNet [139], and others [9, 33, 40].
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.
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.
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., Moon-centered, Moon-fixed) frame,
which may be represented in homogeneous coordinates as
A crater is a planar feature, so let the plane Pi containing crater i be described by the 4 × 1 vector πi, such that
Here, ℙn denotes projective space, which is discussed in greater detail in Ref. [51] and Ref. [43].
Within the plane Pi, 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 Pi with a different origin and different coordinate axis orientation (see Fig. 6). Let the center of the ellipse be located at [xc,yc] 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,
Fig. 6 Visualization of ellipse parameterization in Pi
Substituting this result into Eq. 7 yields the implicit equation
where
Observe that Eq. 9 is the generic equation for a conic, which happens to be an ellipse when B^2^ − 4A**C < 0. A simple calculation with Eqs. 10–12 shows that B^2^ − 4A**C = − 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 Pi be given by
which may be written in homogeneous coordinates as
we observe that Eq. 9 may be rewritten as
where
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 k C (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 (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).
Fig. 7 Visualization of simple conic locus (black) and lines belonging to the conic envelope (red)
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 if x¯ is a point on the conic [51]. For a proper conic, one may therefore compute x¯∝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).
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 Pi from the center of the Moon. The crater’s center point is assumed to lie in the crater plane, πiTp¯Mi(c)=0.
The Moon is very nearly a sphere. Therefore, to an excellent approximation, the crater normal is parallel to the selenographic crater center location pMi(c). If the celestial body were an oblate spheriod, the local “up” direction may be computed through an additional simple calculation [12, 15]. Regardless, assuming a spherical Moon, define a local East-North-Up (ENU) coordinate system, such that
where k^T^ = [001] describes the selenographic description of the lunar pole. With the coordinate axes defined, construct the attitude transformation matrix (passive interpretation of a rotation [146]),
that transforms a vector expressed in the ENU frame of crater i (defined as Ei) to the same vector expressed in the selenograhic frame M.
It is acknowledged that the ENU frame is poorly defined at the poles. In this case, we may replace (25) and (26) with a different convention. Any consistent and well-defined convention for defining the basis vectors of the local horizontal will work.
Since u~i~ is normal to the plane and ρi is the perpendicular distance to the center of the Moon, the vector πi describing the plane Pi may be computed as
Now, consider a 3D point ξEi lying on the rim of crater i as expressed in the crater’s ENU frame Ei. This 3D point may be transformed into selenographic coordinates,
which may be compactly written in homogeneous coordinates
Since the point ξ¯Ei lies in the plane Pi, it is evident that
Thus, defining the columns of TMEi according to
we arrive at a relation between the 2D ellipse point x¯Ei in the crater plane [see Eq. 17] and its corresponding 3D point in the selenographic frame
where k^T^ = [001] and HMi is the 3 × 3 matrix
and where S is the 3 × 2 matrix
Consider a quadratic surface in ℙ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 Pi. 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 ℙ2, points and planes are dual in ℙ3 (a fact that should be evident from Eq. 6). It follows therefore, that one may construct a dual quadric, Qi∗, 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 Qi∗ 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.
Fig. 8 Visualisation of the disk quadric that defines all of planes tangent to Di lying in the plane Pi. A few example tangent planes are illustrated in red
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 ℙ9. The constraint that the rank is ≤ 3 defines a hypersurface in the projective space ℙ9 given by the vanishing of the determinant. The set 𝒵 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 𝒵 (and the hypersurface) is 9 − 1 = 8.
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 Ci is a 3D conic (see Eq. 38 and Table 1), we define its perspective projection into an image as Ai=π(Ci), where π:ℙ3∖{r¯M}→ℙ2 and r¯M is the camera location. Thus, we describe the conic locus for Ai within the image as
where u¯T=[uv1] is the homogeneous coordinate representation of the image pixel coordinate [u, v]. We now introduce a coordinate system to make this projection explicit.
Consider a calibrated camera with selenographic position r~M~ and attitude TCM 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 letting x¯C∝ξC, where x¯C is the projected image plane coordinate of the 3D point ξC. With a digital camera, however, we never observe a point x¯C, but instead we observe the 2D pixel coordinate u^T^ = [u**v] (or u¯T=[uv1]) 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 (u¯) is described by a simple affine transformation,
or, more compactly,
where the camera calibration matrix K is a function of the five calibration dx (and dy) is the ratio of the focal length to pixel pitch in the x (or y) direction, α is the detector skewness, and [up,vp] 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 u¯ in the image according to
where PCM 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 Qi∗ is from Eq. 40 and where A~i~ (which may be computed from Ai∗) 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 HCi to be the 3 × 3 matrix describing the homography between the crater plane Pi and the camera’s image plane
Since HCi 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 si~. For example,
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.
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 SL2-action described below was extensively studied in the 19th century [106]. If a(x, y) = a0x^d^ + a1x^d− 1^y + ⋯ + ady^d^ is a binary form of degree d, then a matrix
acts on it by
where a0′,…,ad′ are polynomial expressions in a0,a1,…,ad,p,q,r,s. For d = 2, the discriminant a12−4a0a2 is invariant under the SL2-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 𝒢 on an arbitrary algebraic variety 𝒱.
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 f1,f2,…,fr such that every other polynomial invariant g is a polynomial expression in the fundamental g=G(f1,f2,…,fr) for some polynomial G(x1,…,xr) in r variables. In the language of commutative algebra, the set of all polynomial functions on the variety 𝒱 form a commutative ring, which we write as ℂ[𝒱] (or ℝ[𝒱] if we work over the real numbers). The subring generated by f1,f2,…,fr is the set of all polynomial expressions in f1,f2,…,fr, which we may write compactly as ℂ[f1,…,fr]. We also define the invariant ring as the set of all invariant polynomials, forming a subring denoted as ℂ[𝒱]𝒢.
Now, observe that f1,f2,…,fr are fundamental invariants if and only if ℂ[𝒱]𝒢=ℂ[f1,f2,…,fr]. Hilbert’s Finiteness Theorem [56, 57] states that there is a finite list of fundamental invariants for groups such as GL~n, SLn~, 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 ℂ(𝒱) be the set of all rational functions on 𝒱, and ℂ(𝒱)𝒢 be the set of all invariant rational functions on 𝒱. We assume that 𝒱 is an irreducible variety so that ℂ(𝒱) has the algebraic structure of a field—specifically, it is the quotient field of the ring ℂ[𝒱]. We now see that ℂ(𝒱)𝒢 is a subfield of ℂ(𝒱). The reader should be warned that the quotient field of ℂ[𝒱]𝒢 can be strictly smaller then ℂ(𝒱)𝒢. 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 f1,f2,…,fr form a system of fundamental rational invariants if every other rational invariant g is a rational expression in f1,f2,…,fr; i.e., of the form G1(f1,…,fr)/G2(f1,…,fr) where G1(x1,…,xr) and G2(x1,…,xr) are polynomials with coefficients in ℂ and G2(f1,…,fr)≠0. There always exists a finite system of fundamental rational invariants f1,f2,…,fr. In fact, it is possible to choose r≤1+dim𝒱 (see Remark 1).
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 𝒱ℝ that parametrizes certain geometric configurations. For example, 𝒱ℝ could be a subvariety of the real projective space ℙℝm defined by homogeneous polynomial equations with real coefficients. The same equations define also a projective variety 𝒱ℂ⊂ℙℂm in complex projective space if we view those equations over the complex numbers. We now assume that 𝒱ℝ is Zariski dense in 𝒱ℂ, which means that any complex polynomial that vanishes on 𝒱ℝ also vanishes on 𝒱ℂ. In this case, we can view the field ℝ(𝒱ℝ) of real rational functions on 𝒱ℝ as a subfield of ℂ(𝒱ℂ). Moreover, the real and complex part of a rational function in ℂ(𝒱ℂ) are rational functions on 𝒱ℝ. This gives a decomposition
where i=−1. We also assume that 𝒢ℂ is a complex algebraic group acting on the variety 𝒱ℂ, such that the subgroup 𝒢ℝ of real group elements acts on 𝒱ℝ and 𝒢ℝ is Zariski dense in 𝒢ℂ. Then, because 𝒢ℝ⊂𝒢ℂ is Zariski dense, a complex rational function is invariant under 𝒢ℂ if and only if it is invariant under the group 𝒢ℝ. Thus, we find that
In particular, if f1,f2,…,fr are real rational invariant functions that generate the field ℝ(𝒱ℝ)𝒢ℝ over ℝ (i.e., ℝ(𝒱ℝ)𝒢ℝ=ℝ(f1,f2,…,fr)), then they also generate ℂ(𝒱ℂ)𝒢ℂ over ℂ (i.e., ℂ(𝒱ℂ)𝒢ℂ=ℂ(f1,f2,…,fr)).
On the other hand, if ℂ(𝒱ℂ)𝒢ℂ=ℂ(f1,f2,…,fr) for some complex rational invariants f1,f2,…,fr, then the real and complex parts of f1,f2,…,fr generate the field ℝ(𝒱ℝ)𝒢ℝ of real 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 f1,f2,…,fr∈L, then K(f1,f2,…,fr) is the set of all rational expressions in f1,…,fr with coefficients in K. This is also called the field that is generated by f1,f2,…,fr over K. So by definition, invariant rational functions f1,f2,…,fr are fundamental rational invariants for the action of 𝒢 on 𝒱 if and only if ℂ(f1,f2,…,fr)=ℂ(𝒱)𝒢.
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 f1,f2,…,fr lie in a field, then we say that f1,f2,…,fr are algebraically dependent over a subfield K if there exists a nonzero polynomial P with coefficients in K with P(f1,f2,…,fr)=0. The transcendence degree of L over the subfield K is the supremum over all r for which there exists f1,f2,…,fr∈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 [62, Theorem VI.1.11]):
If 𝒱 is a variety, then trdeg(ℂ(𝒱)/ℂ) is equal to the dimension of 𝒱. If s=dim𝒱 then we can choose f1,f2,…,fs algebraically independent. If L=ℂ(f1,f2,…,fs) then we have
This shows that ℂ(𝒱)/K is algebraic; i.e., every element of ℂ(𝒱) is algebraic over K. If ℂ(𝒱)≠K, then there exists an element fs+1∈ℂ(𝒱) with
by the Theorem of the Primitive Element (see [72, Theorem 4.6]).
We are particularly interested in trdeg(ℂ(𝒱)𝒢/ℂ) which is the maximal number of rational invariants that are algebraically independent (over ℂ). Let 𝒢⋅v={g⋅v∣g∈𝒢} be the orbit of v∈𝒱 under the 𝒢-action and let 𝒢v={g∈𝒢∣g⋅v=v} be the stabilizer of v. Then we have dim𝒢=dim(𝒢⋅v)+dim𝒢v. If s=maxv∈Vdim(𝒢⋅v) is the largest dimension of an orbit, then almost all orbits in 𝒱 have dimension s and we call s the generic dimension of an orbit in 𝒱. The following results follows from Rosenlicht’s Theorem [122].
If s is the dimension of a generic orbit in 𝒱, then the maximum number of algebraically independent invariants is equal to dim𝒱−s.
Rosenlicht proved in [122] that there exist a 𝒢-stable nonempty (Zariski) open subset 𝒰⊆𝒱 that has a geometric quotient. A geometric quotient is an algebraic variety 𝒰/𝒢 together with a surjective morphism π:𝒰→𝒰/𝒢 such that the fibers of π are exactly all 𝒢-orbits and π also has additional properties such as ℂ(𝒰/𝒢)=ℂ(𝒰)𝒢=ℂ(𝒱)𝒢. The fibers of π have dimension at most s and almost all fibers have dimension s. (Actually one can show that all the fibers have dimension s.) The maximum number of algebraically independent invariants is
□
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 ℙ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 ℙ2. A conic in ℙ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 ℙ5 because such a polynomial has 6 = 5 + 1 coefficients. Let 𝒱 be the variety of d-tuples of conics in ℙ2, Then 𝒱≅(ℙ5)d has dimension 5d. For 𝒢 we take the group PGL3 which has dimension 8. For d = 1 the group PGL3 acts transitively on all nondegenerate quadratic forms. This means that PGL3 has a Zariski dense orbit in 𝒱=ℙ5, so the dimension of a generic orbit is s = 5 and the number of independent rational invariants is dim𝒱−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𝒢=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𝒱−s=5d−8.
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 ℙ3.
Let 𝒵 be the 8-dimensional variety introduced at the end of Section 4.1.2 parameterizing conics in ℙ3, and let 𝒱=𝒵d be the space of d conics in ℙ3. The 15-dimensional group 𝒢=PGL4 acts on ℙ3, 𝒵, and 𝒱. If d = 1 then the group acts transitively on all nondegenerate conics so that 𝒱=𝒵 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 C1 and C2. If the conics are generic enough, then the stabilizer of the pair (C1,C2) is finite. To see this, suppose that g ∈PGL4 lies in the connected component of the identity in the stabilizer of (C1,C2). Then g must also fix the planes P1 and P2 through C1 and C2 respectively. Further, g must also fix the intersections C2∩P1={c¯2(1),c¯2(2)} and C1∩P2={c¯1(1),c¯1(2)} and because g lies in the connected component of the identity, it fixes all the points c¯1(1),c¯1(2),c¯2(1),c¯2(2) on the line P1∩P2 individually. There are two points r¯1(1),r¯1(2) on C1 such that the tangent lines to C1 at r¯1(1) and r¯1(2) go through c¯2(1). So g must also fix r¯1(1) and r¯1(2). Among r¯1(1),r¯1(2),c¯1(1),c¯1(2) in the plane P1 there are no three points on the line (Fig. 9), because they all lie on the same conic C1. This implies that g must fix all the points in the plane P1. Similarly it must fix all the points in the plane P2. From this follows that g is the identity element in PGL4. Thus, the dimension of a general orbit is equal to dimPGL4=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𝒱−15=8d−15.
Fig. 9 Since r¯1(1),r¯1(2),c¯1(1),c¯1(2) all lie on C1 no three of these points may be on the same line
This result agrees with the past literature that has identified one independent rational invariant for a pair of conics in ℙ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.
In Section 5.3.1, we found the invariants for a d-tuple of conics in ℙ3. These, however, are not the invariants we need because we do not measure the crater rims directly in ℙ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.
We begin by discussing the general setup for computing invariants from projections. To do this, we consider models M in ℂ3 space that are parameterized by some variety 𝒱. For example, 𝒱 can consist of all d-tuples of points in 3D, or all d-tuples of conics. We also have camera position o¯∈ℂ3 and a perspective projection π:ℂ3∖{o}→ℙ2. We are interested in rational functions on 𝒱 (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 o¯. We can extend ℂ3 to the projective space ℙ3 and also extend the camera projection to a map ℙ3∖{o¯}→ℙ2, where o¯ is equal to o∈ℂ3⊂ℙ3, but viewed in ℙ3. The group PGL4 acts on ℙ3 and contains the group of affine transformations of ℂ3. The group PGL4 does not act on ℂ3, but it acts on the field of rational functions on ℂ3, because ℂ3 and ℙ3 have the same rational functions. We assume that the class of models parameterized by 𝒱 is closed under the action of PGL4 (for example, PGL4 takes conics to conics). Then PGL4 also acts on the rational functions on 𝒱.
Suppose that f is a rational function on 𝒱 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 o¯ of the camera. So there exists a function h such that f(M)=h(π(M)) for every choice of o¯. Let 𝒢o¯ be the image of the group
in PGL4. An element g∈𝒢o¯ fixes o¯ and all the lines through o¯, so we have
So f is invariant under the action of 𝒢o¯. Since this is true for every choice of o¯, we see that f is invariant under the group 𝒢 generated by all 𝒢o¯, o¯∈ℙ3. One can verify that 𝒢=PGL4. To see this, for example, we observe that the group 𝒢 is closed under conjugation with elements in PGL4. This implies that 𝒢 is a non-trivial normal subgroup of PGL4. Since PGL4 is known to be a simple group, we get 𝒢=PGL4. 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 PGL4.
The simplest example of invariants from projections is the case of a d-tuple of points in ℙ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 ℙ3 but within the context of a simpler (and familiar) example.
Therefore, for example, suppose we have d points in ℙ3. For d ≥ 5 there are 3d − 15 independent invariants for the action of PGL4 (because for d ≥ 5, the generic stabilizer in (ℙ3)d is finite). Perspective projection for some fixed cameragives d points in ℙ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∈(ℙ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 PGL4. The variety of the possible images from the model as we move the camera cuts out a variety 𝒰⊆(ℙ2)d of dimension at most dimPGL4=15. For d ≥ 8, the dimension of 𝒰 is strictly smaller than the dimension of (ℙ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 𝒰 of possible images is not an orbit for any group action on (ℙ2)d: the group PGL4 does not act on (ℙ2)d, and 𝒰 may be bigger than any PGL3 orbit.
To proceed, suppose that 𝒱=(ℙ3)d and 𝒴=(ℙ2)d, o¯∈ℙ3 is the position of the camera, and π:𝒱→𝒴 is the projection. (The map π is only defined for elements in (ℙ3∖{o¯})d.) If h∈ℂ(𝒴) is a rational function, then we can pull it back to get a rational function π⋆h=h∘π∈ℂ(𝒱), where ∘ is the composition. Using the inclusion π⋆:ℂ(𝒴)↪ℂ(𝒱) we may view ℂ(𝒴) as a subfield of ℂ(𝒱). We are interested in rational invariants that can be computed from the image in 𝒴. Such invariants are exactly elements in the intersection field ℂ(𝒴)∩ℂ(𝒱)𝒢. We get the following diagram of field extensions with their transcendence
We will see that ℂ(𝒴)∩ℂ(𝒱)𝒢=ℂ, which means that that there are no non-constant rational invariants that can be computed from the image. Suppose that f∈ℂ(𝒴)∩ℂ(𝒱)𝒢, or to be more precise, there exists an h∈ℂ(𝒴) such that f=h∘π∈ℂ(𝒱)𝒢. We will show that f (and h) must be constant.
To understand what is happening, we introduce an equivalence relation ∼ on (ℙ2)d, the variety of possible images of d points in ℙ2. We say I1∼I2 is true when there exists a model M such that I1,I2 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 (I1∼I2 ⇔ I2∼I1). However, it does not satisfy the transitivity axiom (if I1∼I2 and I2∼I3, then I1∼I3). 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=I0,I1,I2,…,Ir=I′ such that I0∼I1,I1∼I2,…,Ir−1∼Ir. If I≡I′ and I0,I1,…,Ir are as above, then f(I)=f(I0)=f(I1)=⋯=f(Ir)=f(I′).
We will now show that I≡I′ for all images I,I′∈𝒴. Let (p¯1,p¯2,…,p¯d)∈(ℙ2)d=𝒴 and q¯1∈ℙ2. Define o¯∈ℙ3 as the position of a camera and π:ℙ3→ℙ2 as the camera projection. Let L1,L2,…,Ld,N1 be the lines through o¯ in ℙ3 corresponding to the points p¯1,p¯2,…,p¯d,q¯1 respectively. Now, as shown in Fig. 10, define P as the plane through the lines L1 and N1 (which intersect at o¯). Let us choose another camera position o¯′ in the plane P randomly. We randomly choose points a¯1,a¯2,…,a¯d on the lines L1,ℒ2,…,ℒd respectively (i.e., the points a¯1,a¯2,…,a¯d are in general position). Let ℒj′ be the line through o¯′ and a¯j for all j. Since o¯′ and a¯1 lie in the plane P, so does the line L1′. Therefore, the lines L1′ and N1 intersect at some point b¯1. Now we have (p¯1,p¯2,…,p¯d)=π(a¯1,a¯2,…,a¯d) and (p¯1′,p¯2′,…,p¯d′)=π′(a¯1,a¯2,…,a¯d) so (p¯1,p¯2,…,p¯d)∼(p¯1′,p¯2′,…,p¯d′). Moreover, we have (p¯1′,p¯2′…,p¯d′)=π′(b¯1,a¯2,…,a¯d) and (q¯1,p¯2,…,p¯d)=π(b¯1,a¯2,…,a¯d), so (p¯1′,p¯2′,…,p¯d′)∼(q¯1,p¯2,…,p¯d). From this follows that (p¯1,p¯2,…,p¯d)≡(q¯1,p¯2,…,p¯d). Repeating this argument shows that for any points p¯1,p¯2,…,p¯d,q¯1,q¯2,…,q¯d∈ℙ2 we have
So we have (p¯1,p¯2,…,p¯d)≡(q¯1,q¯2,…,q¯d). We conclude that h(p¯1,p¯2…,p¯d)=h(q¯1,q¯2,…,q¯d). In other words, h and f are constant.
Fig. 10 Visualization of geometric quantities used to show the non-existence of invariants from the projection of arbitrary points in ℙ3. Illustrated here are the line-point sets L1∩L1′=a¯1, N1∩L1′=b¯1, L2∩L2′=a¯2, and L3∩L3′=a¯3 that result from the two camera locations o¯¯ and o¯¯′. All quantities in black lie in the plane P
The corollary to a d-tuple of arbitrary points in ℙ3 is a d-tuple of arbitrarily placed conics in ℙ3. We will show that there are no non-constant invariants for d conics in ℙ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 ℙ2. We use an explicit geometric construction to show that all d-tuples are equivalent.
As before, let 𝒵 be the variety of conics in ℙ3 defined at the end of Section 4.1.2. Then, let 𝒱=𝒵d be the variety of d-tuples of conics in ℙ3 and 𝒴=(ℙ5)d be the variety of d-tuples of conics in ℙ2. We define a relation ∼ on 𝒴 by (A1,A2,…,Ad)∼(A1′,A2′,…,Ad′) if there exists two camera projections π and π′ and conics (C1,C2,…,Cd) with π(C1,C2,…,Cd)=(A1,A2,…,Ad) and π′(C1,C2,…,Cd)=(A1′,A2′,…,Ad′). Let ≡ be the equivalence relation generated by ∼.
Suppose that B1,A1,A2,…,Ad⊆ℙ2 are conics. Let ℓ^(1)^,ℓ^(2)^ be two lines in ℙ2 that are tangent to both A1 and B1 with A1,B1 within the same region of ℙ2∖(ℓ(1)∪ℓ(2)) (see Fig. 11).
Fig. 11 Visualization of geometric quantities used to show the non-existence of invariants from the projection of arbitrary conics in ℙ3. The conics and lines in this graphic all lie within the image plane. The backprojection of these planar quantities into ℙ3 form quadric cones (for A1 and B1) or planes (for ℓ^(1)^ and ℓ^(2)^)
We denote the point at which ℓ^(j)^ is tangent to A1 (respectively ℬ1) by a¯1(j) (respectively b¯1(j)). Also, let q¯ be the intersection point of ℓ^(1)^ and ℓ^(2)^. Choose a point o¯∈ℙ3 (the center of the camera) and let π:ℙ3∖{o¯}→ℙ2 be the camera projection. Let Xi=π−1(Ai)∪{o¯} be the cone with top o¯ that corresponds to Ai⊂ℙ2. Also, define Y1=π−1(B1)∪{o¯}. The lines ℓ^(j)^ correspond to the plane L(j)=π−1(ℓ(j))∪{o¯} for j = 1,2. The planes L(1) and ℒ(2) intersect in the line Q=π−1(q¯)∪{o¯}, as shown in Fig. 12. We choose another point o¯′ on the line Q and π′:ℙ3∖{o¯′}→ℙ2 is the camera perspective projection with respect to the camera position o¯′. We choose some random plane P in ℙ3. Then Cj=P∩Xj is a conic in ℙ3 and π(Cj)=Aj. Let Xj′ be the cone with top o¯′ through Cj. By construction, we have π(C1,C2,…,Cd)=(A1,A2,…,Ad) and π′(C1,C2,…,Cd)=(A1′,A2′,…,Ad′) for some conics A1′,…,Ad′ in ℙ2. We will construct a conic D1 in ℙ3 with π(D1)=B1 and π′(D1)=A1′. Note that X1′ is a cone with top o¯′ that is tangent to the planes L(1) and L(2). The intersection S(j)=Y1∩L(j)=π−1(b¯1(j))∪{o¯} is a line through o¯, and T(j)=X1′∩L(j) is a line through o¯′ for j = 1,2. Let r¯(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¯1(j). Moreover we choose a point r¯(3)≠r¯(1),r¯(2) that lies on Y1∩X1′. Let R be the plane through r¯(1),r¯(2),r¯(3). Define D1=R∩X1′. Then D1 is tangent to the planes L(1) and L(2) at the points r¯(1) and r¯(2) respectively. Now π(D1) and B1 both are tangent to ℓ^(1)^ and ℓ^(2)^ at the points b¯(1)=π(r¯(1)) and b¯(2)=π(r¯(2)) respectively, and both contain the point π(r¯(3)). It follows that π(D1)=B1 because a conic in ℙ2 is determined by 3 points and the tangent lines at 2 of the points. One would expect the intersection Y1∩X1′ 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 D1 is one of them.
Fig. 12 Visualization of geometric quantities used to show the non-existence of invariants from the projection of arbitrary conics in ℙ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)
It follows that (B1,A2,A3,…,Ad)=π(D1,C2,…,Cd) and (A1′,A2′,…, Ad′)=π′(D1,C2,…,Cd), and therefore (B1,A2,…,Ad)∼(A1′,A2′,…,Ad′). If we combine this with (A1′,A2,…,Ad′)∼(A1,A2,…,Ad) then we get (B1,A2,…,Ad)≡(A1,A2,…,Ad). Repeating the argument shows that
for all conics A1,A2,…,Ad,B1,B2,…,Bd. This implies that there is no nonconstant rational invariant for d conics in ℙ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).
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 𝒱 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.
First consider the case d = 2. Suppose that C1 and C2 are two conics that lie in the planes given by the linear equations f1 = 0 and f2 = 0 respectively. Then C1 and C2 lie in the degenerate quadric surface f1f2 = 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 + t**f1f2 = 0 contains both conics. If two conics C1 and C2 lie on a nondegenerate surface defined by g = 0, then C1 is defined by f1 = g = 0 and C2 is defined by f2 = g = 0. The intersection C1∩C2 is defined by f1 = f2 = g = 0 and consists of two (possibly complex) points. On the other hand, if C1∩C2 consists of two points, one can show that both conics lie on a nondegenerate quadric surface. To find a pair (C1,C2) in 𝒱, we can choose the conic C1 arbitrarily in the 8-dimensional variety of conics, we choose the plane P2 that contains C2 in the 3-dimensional space of hyperplanes, Now the conic C2 is an element of the 5-dimensional variety of conics in P2, but there are two constraints because the conic has to go through the two points of C1∩P2. So the dimension of the variety 𝒱 is 8 + 3 + (5 − 2) = 14. The 15-dimensional group PGL4 acts on 𝒱. Let H be the connected component of the identity in the stabilizer of (C1,C2)∈𝒱. The conics C1 and C2 intersect in two (possibly complex) points c¯1 and c¯2. So g will fix the points c¯1 and c¯2. If L is the line through c¯1 and c¯2 then g will map L to itself. If c¯3 is a 3rd point on L then g may map c¯3 to another point of L. However, if g also fixes the point c¯3 then a similar argument as in Section 5.3.1 shows that g must fix the planes P1 and P2 pointwise, and must be the identity. This shows that the stabilizer of the pair (C1,C2) 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𝒱=14, this implies that there are no rational invariants for a pair of conics lying on a nondegenerate quadric surface.
Let us now assume that d ≥ 3. If a quadric surface S contains 3 distinct conics C1,C2,C3 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 Pi be the plane through Ci given by fi = 0 where fi is a linear function. Let a¯∈ℙ3. We can multiply h and h′ with nonzero scalars such that (h−h′)(a¯)=0. In the plane Pi, the restriction of the quadratic polynomial h−h′ vanishes on Ci∪{a¯}. This implies that h−h′ is divisible by fi for i = 1,2,3. But then h−h′ is divisible by f1f2f3 and h−h′ must be zero.
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 ℙ3. So the variety 𝒱 of d-tuples of conics that lie on a common quadric surface has dimension 9 + 3d. The stabilizer of a generic point in 𝒱 is finite, so the dimension of a generic orbit is 15=dimPGL4. So the number of independent rational invariants is (9 + 3d) − 15 = 3d − 6 for d ≥ 3.
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.
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 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.
Fig. 13 Two arbitrarily placed non-coplanar conics do not intersect one another. The common 3D line Lij 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 Lij 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
Suppose we have crater i described by the 3D conic Ci that lies in plane Pi. Let the Moon-centered quadric cone of crater i be given by Xi as described in Eq. 37. By construction, this cone must pass through the conic locus Ci, and we interpret Ci as the conic section from Eq. 38.
Now, suppose we have two crater i and crater j. We compute their intersection by substitution of Eq. 38 as,
Define the 3D line Lij as the intersection of planes Pi and Pj,
such that
The intersection of a line with a quadric cone, e.g. Lij∩Xi, produces two points. Define the intersection points of this line with crater i as {ci(1),ci(2)}. 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 {ci(1),ci(2),cj(1),cj(2)} are colinear since they all lie on the line Lij. They are also distinct, meaning their cross ratio ρi**j is a 3D invariant of the conic pair [114]
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 Lij or its intersection points {ci(1),ci(2),cj(1),cj(2)} 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.
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 Pi. 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 Pi with the Moon-centered quadric cone Xi. However, since the crater must lie on both Xi and S, it is clear that we can also write the 3D crater conic as,
Therefore, we see that
Suppose we have two Ci and Cj. Because Ci and Cj are formed by the intersections of planes Pi and Pj with the quadric surface S, we observe that,
where Lij is the 3D line formed by the intersection of the two crater planes (see Eq. 60). This is illustrated in Fig. 14.
Fig. 14 Two non-coplanar conics Ci and Cj 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 Lij formed by the intersection of the two crater planes. The intersection points also lie in the plane Pij formed by the join of Lij and the body center
If the two craters physically intersect, the line Lij pierces S and the intersection points are real. If the craters do not physically intersect, the intersection occurs over the complex numbers. Let s¯ij(1) and s¯ij(2) be the two intersection points
Therefore, define Pij as the plane passing through the center of the Moon and the line Lij (or, equivalently passing through the center of the Moon, s¯ij(1), and s¯ij(2)). Clearly, the numerical stability of computing Pij in this manner is poor if the craters are close to one another and Pi and Pj 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 Xi∩Xj 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 find the other four intersection points are a set of repeated intersections. To see this, combine the above results to find,
The advantage of using the quadric cone intersections to describe s¯ij(1) and s¯ij(2) 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 points s¯ij(1) and s¯ij(2) lie on Ci, the projection of these points lie on the projection of Ci (which we called Ai). The same holds for Cj, and its projection Aj. Therefore,the projection of s¯ij(1) and s¯ij(2) is recoverable as two of the intersection points of the image conics. That is,
where PCM is the projection matrix from Eq. 46 and the pixel coordinate u¯ij(k) is the projection of s¯ij(k) into the image.
The challenge, therefore, is to determine which of the four intersection points of Ai and Aj are the projection of s¯ij(1) and s¯ij(2). There are six possible combinations of these four points, corresponding to six possibilities for the projection of line Lij into the line ℓi**j. We now discuss how to uniquely compute ℓi**j.
Consider the two conics Ai and Aj 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 ℓi**j 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 ℓi**j, which is the projection of Lij. As the intersection of two non-coplanar planes, Lij is real valued by construction. Thus, its projection ℓi**j 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~i**j~ has rank 2. The remaining steps are now to find the lines described by B~i**j, determine if they are real valued, and then determine which one is ℓi**j~.
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~i**j~. 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
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 −Bij∗ are the squares of the elements of z¯. Thus, defining the columns and elements of Bij∗ as
any column of Bij∗ may be scaled to find z¯ according to
where any k ∈{1,2,3} with bk**k≠ 0 will work. Best numerical performance is achieved by selecting the value for k that yields maxk∥bkk∗∥. 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 Bij∗ negative.
Therefore, with z¯ 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 ℓi**j. The other does not. Each of the lines g and h divide ℙ2 into two regions. Since the line we seek comes from the projection of Lij=Pi∩Pj, 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 ℓi**j.
If we were to compute the intersection of ℓi**j with either Ai or Aj we would obtain two points that are the projection s¯ij(1) and s¯ij(2). Fortunately such a computation is not necessary since we work directly with ℓi**j in subsequent discussions.
Briefly, we observe that Ci and the two points s¯ij(1) and s¯ij(2) all lie in the plane Pi. While a projective invariant exists for two points and a conic (all coplanar) [41], the points must not lie on the conic. Since s¯ij(1) and s¯ij(2) lie on Ci by construction, we seek an alternative way of constructing an invariant.
Consider instead the 2D line ℓi**j formed by the projection of 3D line Lij (or, equivalently, by the join of image points u¯ij(1) and u¯ij(2)). 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.
Suppose a projective camera observes three Ci, Cj, and Ck. The intersection of the corresponding planes (Pi, Pj, and Pk) produces the three 3D lines Lij, Lik, and Ljk. 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 aijk, is the apex of the pyramid formed by the three planes,
Now, let the three 3D craters Ci, Cj, and Ck project to 2D conics Ai, Aj, and Ak in the image. For any given pair of 3D conics (e.g., Ci and Cj), the 3D line through their intersection points (e.g., Lij) is coplanar with both of these 3D conics and it is possible to recover its projection (e.g., ℓi**j) from just the projected image conics (e.g., Ai and Aj). Thus, for a triad of observed craters in an image, we may compute the three lines ℓi**j, ℓi**k, and ℓj**k. Since Ci, Lij, Lik are all coplanar, an invariant exists in the image using Ai, ℓi**j, and ℓi**k. The same is true for craters Cj and Ck.
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 ℙ2, we can also form a cross-ratio of four lines passing through a point. Therefore, as shown in Fig. 15, consider a conic Ai. Let the two lines ℓi**j and ℓi**k be described by the 3 × 1 vectors ℓi**j and ℓi**k, and let these two lines intersect at the point p¯=ℓij×ℓik. Using the pole-polar relation for a conic, we may find the two lines from p¯ that are tangent to the conic Ai, which we call w1 and w2. Since these four lines go through the point p¯, we may form a cross ratio that is a projective invariant
Fig. 15 Visualization of contributing components to invariant for a conic Ai and two lines ℓ
i**jand ℓj**k. The lines w1and w2pass through p and are tangent to Ai. The line m is the polar of point p¯ with respect to Ai
We can now form a classical Cayley-Klein metric as
which, from hyperbolic geometry, is known to be equivalent to
where the adjugate matrix Ai∗ describes the conic envelope of Ai. To understand the need for the absolute value in the numerator, recall from Section 4 that ℓi**j and −ℓi**j describe the same line (same for ℓi**k and −ℓi**k). Likewise, Ai∗ and −Ai∗ 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 coshd(ℓij,ℓik)≥1 for all d(ℓi**j,ℓi**k). Absolute values are not necessary in the denominator since this is quadratic in ℓi**j, ℓi**k, and Ai∗. 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 αi2 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(ℓi**j,ℓi**k). The spacecraft navigator will immediately see that this is analogous to how small inter-star 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 (ℓi**j,ℓi**k), 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 Ji, Jj, and Jk are independent. Since we know from Section 5.5 that a d-tuple 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 Ji, Jj, and Jk are independent then we have found all the algebraically independent invariants that exist. This is now shown.
We observe that the invariants αi2,αj2,αk2 (and hence also αi,αj,αk and Ji, Jj, Jk) 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 = t1, y = t2 and z = t3 respectively. Within the plane z = t3 we have the crater given by x2+y2=1−t32 and the lines x = t1 and y = t2. So we get
Thus, we may compute α32
and, by symmetry, we also find
It is possible to verify that the Jacobi matrix
is invertible for some choice of t1,t2,t3 (e.g., t1 = t2 = t3 = 1/2). This implies that α12,α22,α32 (and, therefore, Ji,Jj,Jk) are algebraically independent [75].
Despite the rather long development up to this point, the practical implementation of this method for computing the invariants Ji,Jj,Jk is straightforward and is summarized in Algorithm 1. The reader should note that the invariants Ji,Jj,Jk are computed directly from the measured image conics A1,A2,A3 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 Aj and Ak, 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 J1, J2, and J3 are identical in both images since they are formal projective invariants. This makes evident the power of such invariants for crater pattern recognition.
Fig. 16 Example of invariants J
i, Jj, and Jk(blue) for the same triad of craters (Ai, Aj, and Ak) 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]
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.
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]
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 ℙ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 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.
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
Therefore, let {a¯k}k=14=Ai∩Aj be the four intersection points of the conics in the image plane and let p¯ be an arbitrary point on conic Ai. Since the points are written in homogeneous coordinates, we may form the line joining p¯ and a¯k as wk=p¯×a¯k. Thus, the cross ratio Crw1,w2,w3,w4 for ellipse Ai is a projective invariant, and is the same for any choice of p on the ellipse. A second invariant may be computed in exactly the same way by picking an arbitrary point on ellipse Aj. This cross-ratio may also be rewritten in terms of the eigenvalues of the matrix Ai−1Aj [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 Ci projects to image ellipse Ai 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 si~ and sj. 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, si, or sj~ (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 Ci and Cj. 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| = |Aj~| = 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| = |Aj~| = 1, the two unique invariants are simply
We note that [41] follows a completely different approach to finding Ii**j and Ij**i, 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.
It is established in Section 5.2 that there exist 5d − 8 projective invariants for a d-tuple 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, Aj, and Ak~ such that |A~i| = |Aj| = |Ak| = 1. Thus the coefficients Θi1, Θ4, and Θ7 become trivial, and the remaining seven coefficients become the unique invariants. That is, with |A| = |Aj| = |Ak~| = 1, we find that
Therefore, if three coplanar craters Ci, Cj, and Ck project into three image conics Ai, Aj, and Ak (described by the matrices A~i, Aj, and Ak~), 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.
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 projective and permutation (p^2^) invariants. These conventions each have their advantages, which we now discuss.
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.
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
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 three-element
where the individual elements are from Eqs. 94 to 96.
Likewise, for a triad of conics on a plane, we build a seven-element
where the individual elements are from Eqs. 128 to 131.
As before, suppose we observe a triad of craters. Now, instead of having an order-dependent descriptor, suppose that we want to match to a precomputed index in 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 Jijk′T=[Ji′,Jj′,Jk′], where Ji′≤Jj′≤Jk′ is a rearrangement of (Ji,Jj,Jk), i.e., we sort the numbers Ji,Jj,Jk 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 (Ji′,Jj′,Jk′) is a cyclic permutation of (Ji,Jj,Jk) with Ji′=min(Ji,Jj,Jk).
Thus, for a triad of conics on a quadric surface, we build a three-element
where we label observed crater i such that Ji′=min(Ji,Jj,Jk) 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
where we label observed crater i such that Iij′=min(Iij,Ijk,Iki) 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.
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.
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 Ji,Jj,Jk directly from an image using Eqs. 94 to 96. Our goal is to compute rational functions of Ji,Jj,Jk 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, S3, 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 ℤ/3ℤ, which we may also view as the alternating group A3, we require one more
(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
where F1,F2,F3 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 ℤ/3ℤ act on ℝ3 by permuting the coordinates. Then we desire F1,F2,F3 to be invariant functions on ℝ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 ℤ/3ℤ on ℝ3), we can make this very 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 ζ=e2πi/3=−−1+3i2. Then the generator σ of ℤ/3ℤ acts by the diagonal matrix
We define F1, F2, and F3 as rational functions of a, b, c,
The denominator of F2 and F3 is equal to 12[(x−y)2+(y−z)2+(z−x)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 F2(x, y, z) and F3(x, y, z) converge to 0. Consequently, we can extend F2 and F3 to continuous functions on all of ℝ3 by defining them to be 0 whenever x = y = z.
F(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′)∈ℝ3.
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 F2(x,y,z)=F2(x′,y′,z′)=F3(x,y,z)=F3(x′,y′,z′)=0 then we must have x = y = z and x′=y′=z′ and from F1(x,y,z)=3x=3x′=F1(x′,y′,z′) it follows that (x,y,z)=(x′,y′,z′). Now suppose that F2(x, y, z)^2^ + F3(x, y, z)^2^ (and, therefore, F2(x′,y′,z′)2+F3(x′,y′,z′)2) is nonzero. Then b,b′,c,c′ are all nonzero. Since F2(x,y,z)=F2(x′,y′,z′) and F3(x,y,z)=F3(x′,y′,z′) we get
and
This implies that
Therefore b′=ζsb for some s ∈{0,1,2}. Taking complex conjugation, we get that c′=ζ−sc. We conclude that (a, b, c) and (a′,b′,c′) lie in the same orbit ℤ/3ℤ, which means that (x′,y′,z′) is a permutation of (x, y, z). □
One can also verify that the field of rational invariants ℝ(x,y,z)ℤ/3ℤ is generated by F1,F2,F3.
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, (Ii**j,Ij**i),(Ij**k,Ik**j),(Ik**i,Ii**k), and a 7th invariant for the entire triad Iijk. 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 Iijk is already symmetric and invariant to the ellipse label order.
Simply by applying the results of the previous section, we can find six ℤ/3ℤ-invariant functions from Fℓ(Ii**j,Ij**k,Ik**i) and Fℓ(Ij**i,Ik**j,Ii**k). 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 ℤ/3ℤ acts on (x1,y1,z1) and (x2,y2,z2) by the same cyclic permutation of coordinates, and it acts on b1,b2 with a scalar ζ and on c1,c2 with a scalar ζ^2^. We see, therefore, that b1c2 is invariant. We define G1(x1,y1,z1,x2,y2,z2) and G2(x1,y1,z1,x2,y2,z2) be the real and imaginary part of b1c2, which we compute as
We may isolate the real part to obtain G1
and the imaginary part to obtain G2
Suppose that xℓ,yℓ,zℓ,xℓ′,yℓ′,zℓ′∈ℝ for ℓ = 1,2 such that
Then
is a cyclic permutation of
By Lemma 1 we know that (x1′,y1′,z1′) is a cyclic permutation of (x1,y1,z1) and (x2′,y2′,z2′) is a cyclic permutation of (x2,y2,z2), i.e., b1′=ζsb1 and b2′=ζtb2 for some s, t ∈{0,1,2}. It is not yet clear that this is the same permutation, i.e., whether s = t. From G(x1,y1,z1,x2,y2,z2)=G(x1′,y1′,z1′,x2′,y2′,z2′) follows that b1b¯2=b1c2=b1′c2′=b1′b¯2′. On the other hand, b1b¯2=b1′b¯2′=(ζsb1)(ζ−tb¯2)=ζs−tb1b¯2. So either b1 = 0, b2 = 0, or s = t. If t = s, then we are done. If b1 = 0 then x1=y1=z1=x1′=y1′=z1′ and we are also done because (x2′,y2′,z2′) is a cyclic permutation of (x2,y2,z2). The case b2 = 0 goes similarly. □
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
As was noted above, the descriptiveness of such a scheme is not optimal, because this descriptor remains the same when the values of Ii**j,Ij**k,Ik**i cyclically rotate while the values of Ij**i,Ik**j,Ii**k stay the same. This means that 3 different configurations might have the same descriptor.
A second option is to form a nine-element descriptor
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 (x1,y1,z1,x2,y2,z2) 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 ℝ6 by defining the function to be 0 whenever x1 = y1 = z1 or x2 = y2 = z2.
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=ℂ(x1,y1,z1,x2,y2,z2) generated by Fℓ(x1,y1,z1), ℓ = 1,2,3, F1(x2,y2,z2) and Gℓ(x1,y1,z1,x2,y2,z2), for ℓ = 1,2 is equal to K=ℂ(a1,a2,b13,b1c1,b1c2,c1b2). We observe that L = K(b1) and b1 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ℤ/3ℤ and the degree of the extension L/Lℤ/3ℤ is three. This implies that K=Lℤ/3ℤ. In other words, every ℤ/3ℤ-invariant rational function in Ii**j,Ij**k,Ik**i,Ij**i,Ik**j,Ii**k is a rational function in the first 6 elements of the descriptor.
However, in the degenerate case Ii**j = Ij**k = Ik**i = t for some fixed t we have G(Ii**j,Ij**k,Ik**i,Ij**i,Ik**j,Ii**k) = 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).
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).
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
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.
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 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.3M3=3.66×1017 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].
Fig. 20 Relative lunar crater density as a function of diameter and ellipticity. Areas of darker blue indicate more craters per unit area. Overlayed boxes show regions from which craters are drawn to for the local and global crater indexes
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
First, the lunar surface is tiled into Npix = 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.
Fig. 21 Example lunar crater triad belonging to HEALPix 6318. Overlay is on Apollo Metric Camera image AS17-M-0684. (Credit for raw scans of Apollo flight film NASA/JSC/ASU. See [73, 121]
Fig. 22 Graphical depiction of the number of triads in each HEALPix surface pixel. Left column shows orthographic projection of lunar near side, while the right column shows the same for the lunar far side
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 extensive use for star pattern matching in space applications [95, 134], and the n-d k-vector 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.
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 Npix = 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,7373=1.49×1012 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.
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 non-coplanar 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 Npix = 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 9043=1.23×108 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 Npix = 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 313=4,495 combinations, the HEALPix grouping strategy only produces the 707 crater triads.
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.
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.
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 patterns—thus 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 TCM 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 HCi 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 i th crater’s ENU frame, TCM is presumed known, and k^T^ = [001]. Thus, the only unknown in Eq. 150 is r~M~.
Now, recalling (34) for HMi, 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 rM, and the remaining terms are linear in rM. Therefore, we first find the unknown scalar si~ using only the upper-left 2 × 2 submatrix. Using the Frobenius norm, compute ŝi as
which has the least squares solution
With the scale si 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.
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 (Ai) and the image ellipse we measure of the same crater’s rim (A~i). For the moment, however, we will briefly discuss how to compute the distance between two arbitrary Ai and Aj.
Given two ellipses in an image, Ai and Aj, we seek a scalar distance metric d(Ai,Aj) that satisfies the three usual axioms for a distance metric [32]
as well as a fourth axiom unique to this application
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.
Fig. 23 Examples of an ellipse pair undergoing a simiarity transformation. The left-most ellipse pair may be translated (S
1); translated and rotated (S2); or translated, rotated, and scaled (S3). 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
There exist various ad hoc methods based on the explicit ellipse parameters (a, b, xc, yc, 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(Ai)=1) and then compute the Frobenius norm of their 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 one using the Jaccard distance and another by interpreting the ellipse as a multivariate Gaussian.
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 Ai by finding those pixels that satisfy the inequality u¯TAiu¯<0, which is a simple and parallelizable computation. Once we count the pixels in this way, computing the Jaccard distance is trivial.
It is possible to interpret the 2D ellipse Ai as a bivariate Gaussian probability density function. If that interpretation is used for both Ai and Aj, 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 Ai as a bivariate Gaussian. Therefore, begin by subdividing the conic locus matrix A~i~ describing the image ellipse Ai as
and it follows that that u¯iTAiu¯i=0 and that
With this particular scaling of A~i, we immediately recognize the 2 × 2 submatrix Yi~ as the shape of the ellipse
and yiT=[uci,vci] 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 Ai corresponds to the distribution N(yi,Yi−1) with probability density
Therefore, given two ellipses Ai and Aj, we may compute the distance between these two ellipses as the “angle” between pAi(u) and pAj(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,
where
Thus, noting that c is a constant, we rewrite the integral as
Consequently,
We may therefore compute the distance metric dG**A 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 Ai and Aj.
We propose to use the Gaussian angle distance metric as the criteria for evaluating a potential crater correspondence. In the case where Aj is a perturbed version of Ai, 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.
Fig. 24 Example Monte Carlo (MC) results (50,000 runs) for crater contour correspondence check using Metric Camera images from Apollo 17. The left frame shows results on image AS17-M-0786 collected on 12 December 1972 (17:31:56 UTC), while the right frame shows results on image AS17-M-1557 collected on 13 December 1972 (16:32:57 UTC). The top pair of images show the original Apollo image (left) and with the overlay of a red elliptical crater rim and yellow region of interest (right). The bottom pair of images is the zoomed-in region of interest by itself (left) and with an overlay of the red circular crater rim and 50 yellow example crater rims from the Monte Carlo (right). The perturbed limb fits (e.g, those plotted in yellow) were generated by perturbing the true limb’s ellipse parameters. The histograms at the bottom show the statistic d2/σ2∼χ42. The red overlay on the histogram is the analytic χ42 probability density function (PDF). (Credit for raw scans of Apollo flight film NASA/JSC/ASU. See [73, 121])
Since the distance metric is observed to follow a χ42 distribution, we may construct a statistically-informed criterion for accepting a crater match hypothesis. This is accomplished by evaluating the quantile function for the χ42 distribution. For example, the 99th percentile for the χ42 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 χ42 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 χ42 distribution is modest and numerical experiments can help the analyst decide if this is important for a specific application.
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 (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, uc, vc 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.
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 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.
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.
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.
Fig. 25 Example crater identification on a sampling of images from the Clementine UV/Vis Camera. Each row shows results for a different Clementine image. Left frame is the raw image from [100]. Center and right frame show examples of successful (correct) crater triad matches in white, with centers connected by a yellow triangle to facilitate pattern visualization. Other craters (observed but not contributing to the matched triad) are shown in red. Matches were performed using the seven-element descriptor for local patterns
We demonstrate the three-element global matching descriptor on synthetic images along a reference lunar trajectory.
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, uc, and vc) 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, uc, and vc), 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.
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 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
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 flight-capable system. Crater-based TRN is a rich field for future study.
The authors thank Peter Sturm for insightful discussions related to the projection of conics that led to the technique presented in Section 9.1. The authors also thank Paul McKee and the NASA Johnson Space Center for assistance with production of the PANGU synthetic images. The following individuals provided valuable discussion that greatly improved the quality of this Stuart Robbins, Joseph Mundy, Richard Hartley, Jeffrey Banks, James McCabe, Kristen Schell, Jacob Hinkle, Christopher D’Souza, Juno Woods, and Devin Renshaw. This work was partially supported by NASA award 80NSSC17M0027 (J.A.C.), NSF grant IIS 1837985 (H.D.), NSF grant DMS 2001460 (H.D.), and NASA Lunar Data Analysis Program grant 80NSSC17K0343 (R.W.).
John A. Christian, Email: john.a.christian@gatech.edu
Harm Derksen, Email: ha.derksen@northeastern.edu.
Ryan Watkins, Email: ryan.n.watkins@nasa.gov.